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Abstract 

We  present  a  general  methodology  for  constructing  lattice  Boltzmann  models  of  hydrody¬ 
namics  with  certain  desired  features  of  statistical  physics  and  kinetic  theory.  We  show  how  a 
methodology  of  linear  programming  theory,  known  as  Fourier-Motzkin  elimination,  provides  an 
important  tool  for  visualizing  the  state  space  of  lattice  Boltzmann  algorithms  that  conserve  a 
given  set  of  moments  of  the  distribution  function.  We  show  how  such  models  can  be  endowed 
with  a  Lyapunov  functional,  analogous  to  Boltzmann’s  H,  resulting  in  unconditional  numeri¬ 
cal  stability.  Using  the  Chapman-Enskog  analysis  and  numerical  simulation,  we  demonstrate 
that  such  entropically  stabilized  lattice  Boltzmann  algorithms,  while  fully  explicit  and  perfectly 
conservative,  may  achieve  remarkably  low  values  for  transport  coefficients,  such  as  viscosity. 
Indeed,  the  lowest  such  attainable  values  are  limited  only  by  considerations  of  accuracy,  rather 
than  stability.  The  method  thus  holds  promise  for  high-Reynolds  number  simulations  of  the 
Navier-Stokes  equations. 


Keywords:  computational  fluid  dynamics,  thermodynamics,  hydrodynamics,  entropy,  numerical 
stability,  lattice  gases,  lattice  Boltzmann  equation,  detailed  balance,  complex  fluids 
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Though  lattice  models  have  been  used  to  study  equilibrium  systems  since  the  1920’s,  their  ap¬ 
plication  to  hydrodynamic  systems  is  a  much  more  recent  phenomenon.  Lattice-gas  models  of 
hydrodynamics  were  first  advanced  in  the  late  1980’s,  and  the  related  lattice  Boltzmann  method 
was  developed  in  the  early  1990’s.  This  paper  describes  a  particularly  interesting  and  useful  cate¬ 
gory  of  lattice  Boltzmann  models,  which  we  call  Entropic  Lattice  Boltzmann  models ,  and  provides 
a  number  of  mathematical  tools  for  constructing  and  studying  them.  This  introductory  section  de¬ 
scribes  and  contrasts  lattice  methods  for  hydrodynamics,  traces  their  development,  and  endeavors 
to  place  the  current  work  in  historical  context. 

1.1  Lattice  Gases 

The  first  isotropic  lattice-gas  models  of  hydrodynamics  were  introduced  in  the  late  1980’s  [1,  2,  3]. 
Such  models  consist  of  discrete  particles  moving  and  colliding  on  a  lattice,  conserving  mass  and 
momentum  as  they  do.  If  the  lattice  has  sufficient  symmetry,  it  can  be  shown  that  the  density 
and  hydrodynamic  velocity  of  the  particles  satisfy  the  Navier-Stokes  equations  in  the  appropriate 
scaling  limit. 

This  method  exploits  an  interesting  hypothesis  of  kinetic  theory:  The  Navier-Stokes  equations 
are  the  dynamic  renormalization  group  fixed  point  for  the  hydrodynamic  behavior  of  a  system  of 
particles  whose  collisions  conserve  mass  and  momentum.  We  refer  to  this  assertion  as  a  hypothesis, 
because  it  is  notoriously  difficult  to  demonstrate  in  any  rigorous  fashion.  Nevertheless,  it  is  very 
compelling  because  it  explains  why  a  wide  range  of  real  fluids  with  dramatically  different  molecular 
properties  -  such  as  air,  water,  honey  and  oil  -  can  all  be  described  by  the  Navier-Stokes  equations. 
A  lattice  gas  can  then  be  understood  as  a  “minimalist”  construction  of  such  a  set  of  interacting 
particles.  Viewed  in  this  way,  it  is  perhaps  less  surprising  that  they  satisfy  the  Navier-Stokes 
equations  in  the  scaling  limit. 

For  a  time,  lattice-gas  models  were  actively  investigated  as  an  alternative  methodology  for  com¬ 
putational  fluid  dynamics  (CFD)  [4],  Unlike  all  prior  CFD  methodologies,  they  do  not  begin  with 
the  Navier-Stokes  equations;  rather,  these  equations  are  an  emergent  property  of  the  particulate 
model.  Consequently,  the  use  of  such  models  to  simulate  the  Navier-Stokes  equations  tends  to 
parallel  theoretical  and  experimental  studies  of  natural  fluids.  One  derives  a  Boltzmann  equation 
for  the  lattice-gas  particles,  and  applies  the  Chapman- Enskog  analysis  to  determine  the  form  of  the 
hydrodynamic  equations  and  the  transport  coefficients.  In  fact,  numerical  experimentalists  often 
circumvent  this  theory  by  simply  measuring  lattice-gas  transport  coefficients  in  a  controlled  setting 
in  advance  of  using  them  to  study  a  particular  flow  problem,  just  as  do  laboratory  experimentalists. 

One  often  overlooked  advantage  of  lattice-gas  models  is  their  unconditional  stability.  By  insist¬ 
ing  that  lattice-gas  collisions  obey  a  detailed-balance  condition  x,  we  are  ensured  of  the  validity 
of  Boltzmann’s  H  theorem,  the  fluctuation-dissipation  theorem,  Onsager  reciprocity,  and  a  host 
of  other  critically  important  properties  with  macroscopic  consequences.  By  contrast,  when  the 
microscopic  origins  of  the  Navier-Stokes  equations  are  cavalierly  ignored,  and  they  are  “chopped 
up”  into  finite-difference  schemes,  these  important  properties  can  be  lost.  The  discretized  evolution 
equations  need  no  longer  possess  a  Lyapunov  functional,  and  the  notion  of  underlying  fluctuations 
may  lose  meaning  altogether.  As  the  first  practioneers  of  finite-difference  methods  found  in  the 
1940’s  and  1950’s,  the  result  can  be  the  development  and  growth  of  high-wavenumber  numerical 

lThis  is  certainly  possible  for  an  ideal  single-phase  viscous  fluid.  For  complex  fluids  with  finite-range  interaction 
potentials,  it  is  an  outstanding  problem. 
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instabilities,  and  indeed  these  have  plagued  essentially  all  CFD  methodologies  in  all  of  the  decades 
since.  Such  instabilities  are  entirely  unphysical  because  they  indicate  the  absence  of  a  Lyapunov 
functional,  analogous  to  Boltzmann’s  H.  Numerical  analyists  have  responded  to  this  problem  with 
methods  for  mitigating  these  anomalies  -  such  as  upwind  differencing  -  but  from  a  physical  point 
of  view  it  would  have  been  much  better  if  the  original  discretization  process  had  retained  more  of 
the  underlying  kinetics,  so  that  these  problems  had  not  occurred  in  the  first  place.  Lattice  gases 
represent  an  important  first  step  in  this  direction.  As  was  shown  shortly  after  their  first  appli¬ 
cations  to  hydrodynamics,  lattice-gas  models  for  single-phase  ideal  fluids  can  be  constructed  with 
an  H  theorem  [5]  that  rigorously  precludes  any  kind  of  numerical  instability.  More  glibly  stated, 
lattice  gases  avoid  numerical  instabilities  in  precisely  the  same  way  that  Nature  herself  does  so. 

1.2  Lattice  Boltzmann  Methods 

In  spite  of  these  appealing  features,  the  presence  of  intrinsic  kinetic  fluctuations  makes  lattice-gas 
models  less  than  ideal  as  a  CFD  methodology.  Accurate  values  for  the  velocity  field  at  selected 
locations,  or  even  for  bulk  coefficients  such  as  drag  and  lift,  have  an  intrinsic  statistical  error  that 
can  be  reduced  only  by  extensive  averaging.  On  the  other  hand,  the  presence  of  such  fluctuations 
in  a  simple  hydrodynamic  model  make  lattice  gases  an  ideal  tool  for  those  studying  the  statistical 
physics  of  fluids,  molecular  hydrodynamics,  and  complex-fluid  hydrodynamics;  not  surprisingly, 
these  remain  the  method’s  principal  application  areas  [6].  Consequently,  many  CFD  researchers  who 
appreciated  the  emergent  nature  of  lattice-gas  hydrodynamics  but  wanted  to  eliminate  (or  at  least 
control)  the  level  of  fluctuations,  turned  their  attention  to  the  direct  simulation  of  the  Boltzmann 
equation  for  lattice  gases.  Such  simulations  are  called  lattice  Boltzmann  models  [7,  8,  9,  10,  11]. 
These  models  evolve  the  real-valued  single-particle  distribution  function,  rather  than  the  discrete 
particles  themselves,  and  in  this  manner  eliminate  the  kinetic  fluctuations. 

Early  attempts  along  these  lines  restricted  attention  to  Boltzmann  equations  for  actual  lattice 
gases.  It  was  soon  realized,  however,  that  these  quickly  become  unwieldy  as  the  number  of  possible 
particle  velocities  increases.  In  order  for  lattice  Boltzmann  models  to  become  practical  tools,  it 
was  necessary  to  develop  simplified  collision  operators  that  did  not  necessarily  correspond  to  an 
underlying  lattice-gas  model.  The  most  successful  collision  operators  of  this  type  are  those  of  the 
form  due  to  Bhatnager,  Gross  and  Krook  [12],  and  these  have  given  rise  to  the  so-called  lattice 
BGK  models  [13]. 

1.3  Lattice  BGK  Models 

Lattice  BGK  collision  operators  allow  the  user  to  specify  the  form  of  the  equilibrium  distribution 
function  to  which  the  fluid  should  relax.  For  lattice  gases  obeying  a  detailed  balance  2  condition, 
this  is  known  to  be  a  Fermi-Dirac  distribution.  Having  abandoned  lattice-gas  collision  operators, 
however,  it  seemed  unnecessary  to  continue  to  use  lattice-gas  equilibria,  and  practitioners  exploited 
the  freedom  of  choosing  lattice  BGK  equilibria  to  achieve  certain  desiderata,  such  as  Galilean 
invariance,  and  the  correct  form  of  the  compressible  Navier-Stokes  equations. 

The  alert  reader  will  have  noticed,  however,  that  the  evolution  to  lattice-BGK  methods  has 
jettisoned  the  last  vestiges  of  kinetic-level  physics  that  might  have  been  left  over  from  the  original 
lattice-gas  models.  The  move  to  a  Boltzmann  description  of  the  lattice  gas  eliminated  kinetic 
fluctuations,  but  at  least  it  retained  an  H- theorem;  more  precisely,  its  global  equilibrium  still 
extremized  a  Lyapunov  functional  of  the  dynamics.  The  move  to  lattice-BGK  operators  and  the 

2Actually,  a  weaker  condition  called  semi- detailed  balance  suffices  for  this  purpose. 
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arbitrary  legislation  of  the  equilibrium  distribution  function,  however,  completely  abandoned  even 
the  concept  of  detailed  balance  and  the  H- theorem.  Without  a  Lyapunov  functional,  lattice-BGK 
methods  became  susceptible  to  a  wide  variety  of  numerical  instabilities  which  are  ill-understood 
and  remain  the  principal  obstacle  to  the  wider  application  of  the  technique  to  the  present  day. 

1.4  Motivation  for  Entropic  Lattice  Boltzmann  Methods 

In  this  paper,  we  argue  that  the  most  natural  way  to  eliminate  numerical  instabilities  from  lattice 
Boltzmann  models  is  to  return  to  the  method  its  kinetic  underpinnings,  including  the  notion  of  a 
Lyapunov  functional.  We  present  a  general  program  for  the  construction  of  “entropically  stabilized” 
lattice  Boltzmann  models,  and  illustrate  its  application  to  several  example  problems. 

The  crux  of  the  idea  is  to  encourage  the  model  builder  to  specify  an  appropriate  Lyapunov  func¬ 
tional  for  the  model,  rather  than  try  to  blindly  dictate  an  equilibrium  state.  Of  course,  specifying 
a  Lyapunov  functional  determines  the  equilibrium  distribution,  but  it  also  governs  the  approach  to 
this  equilibrium.  It  can  therefore  be  used  to  control  the  stability  properties  of  the  model.  It  should 
be  emphasized  that  the  presence  of  a  Lyapunov  functional  guarantees  the  nonlinear  stability  of  the 
model,  which  is  a  much  stronger  condition  than  linear  stability. 

As  we  shall  show,  collision  operators  that  are  constructed  in  this  way  are  similar  in  form  to  the 
lattice  BGK  collision  operators,  except  that  their  relaxation  parameter  may  depend  on  the  current 
state.  As  a  result,  the  transport  coefficients  may  have  a  certain  minimum  value  in  models  of  this 
type  but,  as  we  shall  see,  this  value  may  actually  be  zero.  In  fact,  the  ultimate  limitation  to  the 
application  of  these  algorithms  for  very  small  transport  coefficients  will  come  from  considerations 
of  accuracy,  rather  than  stability. 

1.5  Structure  of  Paper 

The  plan  of  this  paper  is  as  follows:  In  Section  2  we  present  the  general  program  of  construction  of 
entropic  lattice  Boltzmann  models.  This  includes  the  introduction  of  a  conservation  representation 
of  the  lattice  Boltzmann  distribution  function,  by  means  of  which  the  collision  process  is  most 
naturally  described.  This  representation  also  provides  an  interesting  geometric  interpretation,  as  it 
allows  us  to  describe  a  collision  process  as  a  mathematical  map  from  a  certain  polytope  into  (and 
perhaps  onto)  itself.  We  then  show  a  way  of  constructing  Lyapunov  functionals  on  this  polytope, 
contrast  these  with  Boltzmann’s  H,  and  use  them  to  construct  collision  operators  for  absolutely 
stable  models. 

Section  3  carries  out  this  program  for  a  very  simple  lattice  Boltzmann  model  of  diffusion  in  one 
dimension.  The  model  is  useful  from  a  pedagogical  point  of  view,  since  an  entropically  stabilized 
collision  operator  can  be  written  for  it  in  closed  form.  Moreover,  it  is  possible  to  determine  the 
transport  coefficient  of  this  model  analytically.  The  result  is  a  fully  explicit,  perfectly  conservative, 
absolutely  stable  method  of  integrating  the  diffusion  equation.  We  present  numerical  simulations 
for  various  values  of  the  diffusivity  in  order  to  probe  the  limitations  of  the  technique.  In  the  course 
of  presenting  this  model,  we  show  that  a  linear  algebraic  procedure,  called  the  Fourier-Motzkin 
algorithm,  is  very  useful  for  constructing  and  visualizing  the  conservation  representation. 

Section  4  applies  the  method  to  a  simple  five- velocity  model  of  fluid  dynamics  in  one  dimension, 
first  considered  by  Renda  et  al.  in  1997  [14].  Here  the  geometric  picture  involves  a  four  dimensional 
polytope,  and  the  Fourier-Motzkin  algorithm  is  shown  to  be  very  useful  in  describing  it.  Indeed, 
Appendix  A  gives  the  reader  a  “tour”  of  this  polytope. 
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Finally,  Section  5  applies  the  method  to  the  two-dimensional  FHP  model  and  the  three- 
dimensional  FCHC  model  of  fluid  dynamics.  For  these  examples,  we  are  able  to  provide  appropriate 
conservation  representations.  Though  the  latter  example  has  too  many  degrees  of  freedom  to  allow 
geometric  visualization  of  the  collision  dynamics,  it  is  nevertheless  possible  to  develop  and  present 
a  Lyapunov  functional  and  an  entropic  model.  We  use  these  more  complicated  models  to  demon¬ 
strate  that  entropic  lattice  Boltzmann  models  are  not  computationally  prohibitive,  even  when  the 
number  of  velocities  involved  is  large.  We  describe  the  computational  aspects  of  such  models  in 
some  detail,  and  end  with  a  discussion  of  Galilean  invariance. 
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2  Entropic  Lattice  Boltzmann  Models 

In  this  section,  we  present  the  general  program  of  construction  of  entropic  lattice  Boltzmann  models. 
We  note  that  the  lattice  Boltzmann  distribution  function  admits  a  variety  of  representations,  and  we 
show  how  to  transform  between  these.  In  particular,  we  introduce  a  conservation  representation ,  by 
means  of  which  the  collision  process  has  a  very  natural  description.  We  also  examine  the  geometric 
aspects  of  this  representation,  and  use  these  insights  to  show  how  to  construct  lattice  Boltzmann 
collision  operators  possessing  a  Lyapunov  functional. 

2.1  Representation 

Lattice  Boltzmann  models  are  constructed  on  a  D  dimensional  regular  lattice  £,  and  evolve  in 
discrete  time  intervals  At.  We  imagine  that  a  population  of  particles  exists  at  each  site  x  6  jC,  and 
that  these  particles  can  have  one  of  a  finite  number  n  of  velocities,  Cj/At  where  j  €  {1, . . . ,  n}.  The 
displacement  vectors  cj  are  required  to  be  combinations  of  lattice  vectors  with  integer  coefficients, 
so  that  if  x  6  £,  then  x  +  Cj  e  C.  We  note  that  there  is  nothing  precluding  some  of  the  Cj  from 
being  the  same,  or  from  being  zero;  indeed,  this  latter  choice  is  often  made  to  incorporate  so-called 
“rest  particles”  into  the  model. 

Mathematically,  the  state  of  the  system  of  particles  is  represented  by  n  real  numbers  at  each 
site  x  e  C  and  time  step  t.  These  are  denoted  by  Nj(x,t)  6  5?,  and  can  be  thought  of  as  the 
expected  number  of  particles  with  mass  mj  and  velocity  Cj/At  at  site  x  and  time  t.  As  such, 
we  expect  them  to  be  nonnegative,  Nj(x,t)  >  0.  Taken  together,  these  quantities  constitute  the 
single-particle  distribution  function  of  the  system.  This  is  all  the  information  that  is  retained  in  a 
Boltzmann  description  of  interacting  particles. 

In  what  follows,  we  shall  think  of  the  set  of  single-particle  distribution  functions  as  a  vector 
space.  That  is,  we  shall  regard  the  n  quantities  Nj(x,t),  for  i  6  {1, . .  .,n},  as  the  components  of 
an  n-vector  or  “ket,”  |N(x,t))  e  5?n.  If  no  ambiguity  is  likely  to  result,  we  shall  often  omit  the 
explicit  dependence  on  t  and  write  simply  |N(x)).  If  discussion  is  restricted  to  a  single  lattice  site, 
we  may  further  abbreviate  this  as  |N). 

In  one  discrete  time  step  At,  the  state  of  the  system  is  modified  in  a  manner  that  is  intended  to 
model  collisions  between  the  particles  at  each  site,  followed  by  propagation  of  the  collided  particles 
to  new  sites.  The  collisions  are  modelled  by  modifications  to  the  distribution  function 

|N(x))  -+  |N')  .  (1) 

Here  and  henceforth,  we  use  a  prime  to  denote  the  postcollision  state.  It  is  usually  the  case  that 
|N'(x))  depends  only  on  |N(x)),  though  in  more  sophisticated  lattice  Boltzmann  models  it  may  also 
depend  on  the  |N)’s  in  a  local  neighborhood  of  sites  about  x.  Following  convention,  we  define  the 
collision  operator  Cy  (|N(x)))  as  the  difference  between  the  old  and  new  values  of  the  single-particle 
distribution, 

Cy(|N(x)»  =  iV'(x)  -  Nj(x).  (2) 

As  this  is  the  difference  of  two  kets,  it  can  also  be  thought  of  as  a  ket, 

|C)  =  |N')-|N>.  (3) 
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2.2  Conservation  Laws 

The  collision  process  specified  in  Eq.  (1)  is  usually  required  to  conserve  some  number  of  locally 
defined  quantities.  Usually,  these  quantities  are  additive  over  the  lattice  sites  and  directions.  For 
example,  we  may  require  that  there  be  a  conserved  mass 

M  =  ( 4 ) 

xe£ 


where  we  have  defined  the  mass  density 

n  n 

p(*.  *)  =  l)  =  Yl  *)•  (5) 

j= 0  j= o 

To  streamline  the  notation  for  this,  we  can  define  a  covector  or  “bra”  by  the  prescription 


(p\j  =  rrij, 


(6) 


and  write  simply 


p={p|N>  =  (p|N/). 


(7) 


Thus,  to  the  extent  that  we  think  of  the  single-particle  distribution  function  as  a  vector,  to  each 
conserved  quantity  there  corresponds  a  covector  such  that  the  value  of  the  conserved  quantity  is 
the  contraction  of  the  two.  Because  collisions  are  required  to  preserve  this  value,  the  contraction 
of  the  covector  with  the  collision  operator  must  vanish, 


(P|C>  =  0, 


(8) 


as  can  be  seen  from  Eqs.  (3)  and  (7). 

Lattice  Boltzmann  models  may  also  conserve  momentum, 

^(x,  t)  =  i)>  =  l]mj^Aj(x,i).  (9) 

j—0  j= 0 

This  is  still  linear  in  the  Nj,  so  that  we  can  define 

M  j  = 

and  write 

tt(x,  t)  =  (tt  |  N)  =  (tt  |  N')  .  (11) 

Here  we  must  be  a  bit  careful  to  define  the  meaning  of  the  quantities  involved:  Note  that  is 
a  covector  in  its  j  index,  but  a  vector  in  its  spatial  index  3.  Thus,  the  contraction  with  |N);-  is 
over  the  j  index,  and  results  in  the  vector  7r.  Once  again,  this  covector  annihilates  the  collision 
operator, 

(*■  I  C)  =  0,  (12) 

where  the  right-hand  side  is  a  null  vector. 

3The  pedant  will  note  that  momentum  is  more  properly  thought  of  as  a  covector  in  its  spatial  index,  but  we  do 
not  bother  to  distinguish  between  spatial  vectors  and  covectors  in  this  paper. 
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The  subsequent  propagation  process  is  described  mathematically  by  the  prescription 


Nj(x  +  cj,  t  +  At)  =  Nj(x,  t).  (13) 

Note  that  this  operation  must  be  carried  out  simultaneously  over  the  entire  lattice.  This  may 
alter  the  values  of  the  conserved  quantities  at  each  site,  but  because  it  is  nothing  more  than  a 
permutation  of  the  values  of  the  distribution  function  about  the  lattice,  it  clearly  leaves  unaltered 
the  global  values  of  the  conserved  quantities, 

.  ExerPOM)  and  Exc£7r(x><)-  - .  (14) 

Combining  Eqs.  (2)  and  (13),  we  find  the  general  dynamical  equation  for  a  lattice  Boltzmann  model, 

Nj(x  +  Cj,  t  +  At)-  Nj(x)  =  C,(|N(x))).  (15) 

For  compressible  fluids,  it  is  also  necessary  to  pay  attention  to  conservation  of  energy.  Conser¬ 
vation  of  kinetic  energy  can  be  expressed  using  the  bra, 

=  (16) 

Indeed,  the  right-hand  side  could  be  generalized  without  difficulty  to  anything  that  depends  only 
on  j.  The  problem  of  including  a  potential  energy  function  between  particles  at  different  sites, 
however,  is  much  more  difficult.  If  we  suppose  that  there  is  a  potential  I^(|x_,-  —  x*|)  between 
particles  at  xy,Xfc  e  £,  then  two  problems  arise:  First,  the  total  potential  energy  will  depend 
on  the  pair  distribution  function  -  that  of  finding  two  particles  a  certain  distance  apart.  This  is 
outside  the  scope  of  Boltzmann  methods,  which  retain  only  the  single-particle  distribution.  We 
can  avoid  this  problem  by  making  the  mean-field,  approximation,  in  which  the  probability  of  having 
one  particle  at  Xj  and  another  at  x*  is  simply  the  product  of  the  two  single-particle  probabilities, 
but  then  the  potential  energy  is  quadratic,  rather  than  linear,  in  the  Nj's.  Second,  and  perhaps 
more  distressingly,  the  potential  energy  is  not  preserved  by  the  propagation  step  [15].  These 
considerations  make  it  very  difficult  to  add  potential  interaction  to  lattice  Boltzmann  models.  If 
we  are  willing  to  give  up  the  idea  of  an  exactly  conserved  energy,  and  instead  consider  isothermal 
systems,  then  methods  of  skirting  these  difficulties  have  been  known  and  actively  investigated  for 
some  time  now  [16,  17].  More  recently,  a  method  of  incorporating  interaction  that  maintains  exact 
(kinetic  plus  potential)  energy  conservation  has  also  been  proposed  [18].  In  any  case,  this  paper 
will  be  restricted  to  systems  with  kinetic  energy  only.  The  possibility  of  extending  our  methods  to 
models  with  nontrivial  interaction  potentials  is  discussed  briefly  in  the  Conclusions  section. 

2.3  Geometric  Viewpoint 

The  requirements  of  maintaining  the  conserved  quantities  and  the  nonnegativity  of  the  distribution 
function  place  very  important  constraints  on  the  collision  process.  Much  of  this  paper  will  be  de¬ 
voted  to  describing  -  algebraically  and  geometrically  -  the  most  general  set  of  collisional  alterations 
of  the  Nj(x)  that  meet  these  requirements. 

Suppose  that  there  are  nc  <  n  independent  conserved  quantities.  As  described  above,  these 
define  nc  linearly  independent  covectors  or  bras,  whose  contraction  with  the  collision  operator 
vanishes.  We  shall  sometimes  adopt  special  names  for  these;  for  example,  that  corresponding  to 
mass  conservation  shall  be  called  (p|,  that  for  momentum  conservation  shall  be  called  (rr|,  etc. 
Generically,  however,  let  us  refer  to  them  as  (ACT|,  where  a  =  1 , . .  ,,nc.  These  covectors  are  not 
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necessarily  orthogonal  in  any  sense.  Note,  for  example,  that  the  covectors  for  mass  and  kinetic 
energy  in  Eqs.  (6)  and  (16)  are  not  orthogonal  with  respect  to  the  Euclidean  metric.  In  fact,  there 
is  no  reason  to  insist  on  any  kind  of  metric  in  this  covector  space. 

The  nc  covectors  are  not  uniquely  defined,  since  a  linear  combination  of  two  conserved  quantities 
is  also  conserved.  The  nc  dimensional  subspace  defined  by  the  covectors,  on  the  other  hand,  is 
uniquely  defined,  and  we  shall  call  it  the  hydrodynamic  subspace4.  It  is  always  possible  to  construct 
n  —  nc  more  covectors  that  are  linearly  independent  of  each  other,  and  of  the  nc  corresponding  to 
the  conserved  quantities.  For  example,  this  may  be  done  by  the  Gram-Schmidt  procedure.  Let  us 
generically  call  these  (a^|,  where  rj  =  1, . . .,  n  —  nc;  if  there  are  only  a  few  of  these  and  we  want 
to  avoid  excessive  subscripting,  we  shall  use  the  successive  Greek  letters  (a| ,  (/3| ,  (7I , . . .  for  this 
purpose.  Once  again,  these  are  not  uniquely  defined  or  orthogonal  in  any  sense,  but  they  do  span 
an  n  —  nc  dimensional  subspace  which  we  shall  call  the  kinetic  subspace.  The  union  of  all  n  of 
our  covectors,  namely  the  (A^l’s  and  the  (a^l’s,  thus  constitute  a  basis  for  the  full  n  dimensional 
covector  space. 

The  ket  |N)  whose  components  are  the  single-particle  distribution  function  is  thus  seen  to  be 
but  one  representation  of  the  dynamical  variable.  As  has  been  recognized  for  some  time  now  [19], 
we  are  free  to  make  a  change  of  basis  in  the  lattice  Boltzmann  equation.  In  fact,  one  other  basis 
suggests  itself  naturally.  Given  a  basis  of  covectors  or  bras,  it  is  always  possible  to  construct  a 
dual  basis  of  vectors  or  kets.  Indeed,  our  decomposition  of  the  covector  space  into  hydrodynamic 
and  kinetic  subspaces  is  naturally  mirrored  by  a  corresponding  decomposition  of  the  vector  space. 
Thus  we  construct  |Aff)  where  a  =  1, . . .,  nc  and  |a^)  where  (3  —  1, . . . ,  n  -  nc,  such  that 

(Ac-  I  Ar)  =  SaT  (Act  I  OCfj)  =  0  (17) 

(cfcjj  |  Act)  —  0  (oijj  |  a$)  =  Sq g. 

Thus,  we  can  expend  |N)  in  terms  of  these  hydrodynamic  and  kinetic  basis  kets. 

nc  n— nc 

|N)  =  ]TAct|Act>+  £o,K>.  (18) 

I  <7=1  »?=1 

We  shall  call  the  set  of  coefficients  \a  and  an  the  conservation  representation.  The  explicit  con¬ 

struction  of  this  representation  is  best  illustrated  by  example,  and  we  give  several  of  these  in  the 
following  sections. 

The  advantage  of  the  conservation  representation  for  describing  collisions  is  clear:  When  |N)  is 
expanded  in  terms  of  hydrodynamic  and  kinetic  kets,  collisions  may  change  only  the  coefficients  of 
the  latter.  From  Eqs.  (17)  and  (18),  we  see  that  the  coefficients  of  the  hydrodynamic  kets,  namely 
the  Act’s,  are  the  conserved  quantities  themselves.  These  are  precisely  what  must  remain  unchanged 
by  a  collision.  So,  in  the  conservation  representation,  the  collision  process  is  simply  an  alteration 
of  the  n  —  nc  coefficients  of  the  kinetic  kets,  namely  the  o^’s.  Thus,  this  representation  effectively 
reduces  the  dimensionality  of  the  space  needed  to  describe  the  collision  process  to  n-nc.  Again,  we 
shall  illustrate  the  construction  of  and  transformations  between  these  alternative  representations  of 
the  single-particle  distribution  function  for  several  examples  of  entropic  lattice  Boltzmann  models. 

4Note  that  we  are  using  the  term  “hydrodynamic”  here  to  describe  degrees  of  freedom  that  will  result  in  macro¬ 
scopic  equations  of  motion,  whether  or  not  they  are  those  of  a  fluid.  If  mass  is  the  only  conserved  quantity,  a  diffusion 
equation  generically  results;  if  momentum  and/or  energy  are  conserved  as  well,  a  set  of  fluid  equations  generically 
results.  We  use  the  term  “hydrodynamic”  in  either  case.  This  terminology  is  standard  in  kinetic  theory  and  statistical 
physics,  but  often  seems  strange  to  hydrodynamicists. 
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2.4  Nonnegativity 

While  the  conservation  representation  of  |N)  is  most  natural  for  describing  collisions,  the  original 
representation  is  more  natural  for  describing  the  constraint  of  nonnegativity  of  the  distribution 
function  components.  In  the  original  representation,  we  had  a  restriction  to  Nj  >  0  for  j  =  1, . . . ,  n. 
In  the  conservation  representation,  these  n  inequalities  transform  to  a  corresponding  set  of  linear 
inequalities  on  the  hydrodynamic  and  kinetic  parameters,  \a  and  ar/.  As  is  well  known,  such  a  set 
of  linear  inequalities  define  a  convex  polytope  in  the  parameter  space.  This  construction,  in  the  full 
n  dimensional  space  of  hydrodynamic  and  kinetic  parameters,  shall  be  called  the  master  polytope. 
Specification  of  the  hydrodynamic  parameters  then  defines  the  cross  section  of  the  master  polytope 
that  bounds  the  kinetic  parameters,  o^;  we  shall  call  these  cross  sections  the  kinetic  polytopes ,  and 
it  is  clear  that  they  must  also  be  convex. 

We  shall  construct  the  master  and  kinetic  polytopes  for  simple  entropic  lattice  Boltzmann 
models  later  in  this  paper.  We  shall  see  that  they  become  very  difficult  to  visualize  when  the 
number  of  particle  velocities  becomes  large.  This  difficulty  raises  the  question  of  whether  or  not 
there  is  a  general  method  to  describe  the  shape  of  polytopes  defined  by  linear  equalities  in  this 
way.  It  turns  out  that  such  a  method  is  well  known  in  linear  programming  and  optimization  theory, 
and  is  called  Fourier- Motzkin  elimination  [20] .  It  is  constructive  in  nature,  and  works  for  any  set 
of  inequalities  in  any  number  of  unknowns.  If  the  inequalities  cannot  be  simultaneously  satisfied, 
the  method  will  indicate  that.  Otherwise,  it  will  yield  an  ordered  sequence  of  inequalities  for  each 
variable,  the  bounds  of  which  depend  only  on  the  previously  bounded  variables  in  the  sequence.  If 
it  were  desired  to  perform  a  multiple  integral  over  the  polytopic  region,  for  example,  the  Fourier- 
Motzkin  method  would  provide  a  systematic  procedure  for  setting  up  the  limits  of  integration. 

2.5  Collisions 

Once  we  are  able  to  construct  and  visualize  these  polytopes,  it  is  straightforward  to  describe  the 
constraints  that  conservation  imposes  on  the  collision  process:  Because  the  propagation  process  is 
nothing  more  than  a  permutation  of  the  values  of  the  Nj(x)  on  the  lattice,  it  is  clear  that  it  will 
maintain  nonnegativity.  That  is,  if  the  iVj(x)  were  positive  prior  to  the  propagation,  they  will  be 
so  afterward.  This  means  that  the  post-propagation  state  of  each  site  will  lie  within  the  allowed 
master  polytope.  Given  such  a  post-propagation  state  |N),  we  transform  to  the  conservation 
represenation.  The  coefficients  Xa  of  the  hydrodynamic  kets  are  the  conserved  quantities  and  must 
remain  unchanged  by  the  collision.  Geometrically,  these  determine  a  cross  section  of  the  master 
polytope  within  which  the  state  must  remain.  This  cross  section  is  the  kinetic  polytope  of  the  pre¬ 
collision  state.  The  essential  point  is  that  the  post-collision  state  must  also  lie  within  this  kinetic 
polytope  to  preserve  nonnegativity.  To  the  extent  that  the  postcollision  state  is  determined  by  the 
precollision  state,  this  means  that: 

Collision  Property  1  The  collision  process  is  a  map  from  the  kinetic  polytope  into  itself. 

We  note  that  this  requirement  is  not  without  some  controversy.  It  may  be  argued  that  lat¬ 
tice  Boltzmann  algorithms  are  ficticious  kinetic  models  from  which  realistic  hydrodynamics  are 
emergent.  Since  the  details  of  the  kinetics  are  ficticious  anyway,  why  not  also  dispense  with  the 
requirement  that  the  single-particle  distribution  function  be  positive?  As  long  as  the  conserved 
densities  of  positive-definite  quantities,  such  as  mass  and  kinetic  energy,  are  positive,  why  should 
we  care  if  the  underlying  lattice  Boltzmann  distribution  function  is  likewise?  There  are  two  rea¬ 
sons:  The  first  is  that,  even  if  the  system  is  initialized  with  nonnegative  physical  densities,  the 
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propagation  step  may  give  rise  to  negative  physical  densities  if  negative  values  of  the  distribution 
are  allowed.  To  see  this,  imagine  a  postcollision  state  in  which  all  of  the  neighbors  of  site  x  have 
a  single  negative  distribution  function  component  in  the  direction  heading  toward  x.  Even  if  all 
these  neighboring  sites  had  positive  mass  density,  site  x  will  have  a  negative  mass  density  after  one 
propagation  step. 

The  second  reason  for  demanding  nonnegativity  of  the  distribution  function  is,  to  some  extent,  a 
matter  of  taste.  We  like  to  think  that  the  kinetic  underpinnings  of  the  lattice  Boltzmann  algorithm 
are  more  than  just  a  mathematical  trick  to  yield  a  desired  set  of  hydrodynamic  equations.  Though 
there  may  be  no  physical  system  with  such  a  collision  operator  (not  to  mention  dynamics  constrained 
to  a  lattice),  we  feel  that  the  more  properties  of  real  kinetics  that  can  be  maintained,  the  more 
useful  the  algorithm  is  likely  to  be.  This  is  particularly  important  for  complex  fluids,  for  which 
the  form  of  the  hydrodynamic  equations  is  often  unknown,  and  for  which  we  must  often  appeal  to 
some  kinetic  level  of  description.  Nevertheless,  we  shall  revisit  this  question  in  Subsection  5.5. 

Collision  Property  1  is  simple  in  statement  and  motivation,  but  in  fact  it  weeds  out  many  puta¬ 
tive  lattice  Boltzmann  collision  operators,  including  those  most  commonly  used  in  computational 
fluid  dynamics  research  today  -  namely,  the  overrelaxed  lattice  BGK  operators.  For  sufficiently 
large  overrelaxation  parameter  (collision  frequency),  such  operators  are  well  known  to  give  rise  to 
negative  values  of  the  Nj(x).  This  is  usually  symptomatic  of  the  onset  of  a  numerical  instabil¬ 
ity  in  the  lattice  Boltzmann  algorithm.  We  shall  discuss  such  instabilities  in  more  detail  later  in 
this  paper.  For  the  present,  we  emphasize  that  we  are  not  claiming  that  Collision  Property  1  will 
eliminate  such  instabilities.  The  property  does  mandate  a  lower  bound  of  zero  on  the  Nj(x),  and 
hence  it  restricts  the  manner  in  which  such  instabilities  might  grow  and  saturate.  Nevertheless,  it 
is  generally  still  possible  for  collision  operators  that  obey  Collision  Property  1  to  exhibit  numerical 
instability  if  they  lack  a  Lyapunov  functional. 

2.6  Reversibility 

Collision  Property  1  is  the  minimum  requirement  that  we  impose  on  our  lattice  Boltzmann  collision 
operators.  It  is  possible  to  define  more  stringent  requirements  for  them,  and  we  shall  continue  to 
do  exactly  that  to  satisfy  various  desiderata.  For  example,  we  might  want  to  demand  that  our 
lattice  Boltzmann  algorithm  be  reversible.  A  reversible  algorithm  could  be  run  backwards  in  time 
from  any  final  condition  by  alternately  applying  the  inverse  propagation  operator, 


Nj(x  -  Cj)  <-  Nj(x),  (19) 

followed  by  an  inverse  collision  operator  to  recover  an  initial  condition  at  an  earlier  time.  For 
such  an  inverse  collision  operator  to  exist,  the  map  from  the  allowed  polytope  to  itself  must  be 
one-to-one.  Since  we  have  already  demanded  that  the  map  be  into,  this  means  that  it  must  also 
be  onto,  hence: 

Collision  Property  2  A  collision  process  is  reversible  if  it  is  a  map  from  the  allowed  polytope 
onto  itself. 

2.7  Imposition  of  a  Lyapunov  Functional 

The  criteria  that  we  have  set  out  thusfar  assure  that  the  conservation  laws  -  and  hence  the  First 
Law  of  Thermodynamics  -  will  be  obeyed.  To  ensure  stability  and  thermodynamic  consistency, 
however,  it  is  necessary  to  also  incorporate  the  requirements  of  the  Second  Law. 
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We  suppose  that  our  system  has  a  function  H  of  the  state  variables  Nj  that  is  additive  over 
sites  x, 

ff=][>(N(x)),  (20) 

x€jC 

and  additive  over  directions  j, 

h(N(x))  =  '£ej(Nj(x)),  (21) 

j=i 

where  the  functions  6j  are  defined  on  the  nonnegative  real  numbers.  One  might  wonder  if  Eq.  (21) 
could  be  generalized,  but  in  fact  it  has  been  shown  to  be  a  necessary  condition  for  an  ^-theorem  [21]. 
It  is  clear  from  this  construction  that  H  is  preserved  under  the  propagation  operation  of  Eq.  (13). 
If  we  require  that  our  collisions  never  decrease  the  contribution  to  H  from  each  site  -  that  is,  that 
h  (N(x))  can  only  be  increased  by  a  collision  -  then  H  is  a  Lyapunov  functional  for  our  system, 
and  the  existence  and  stability  of  an  equilibrium  state  is  guaranteed. 

We  note  that  h  is  a  function  of  the  Nj  for  all  Nj  >  0,  and  can  therefore  be  thought  of  as  a  scalar 
function  on  the  master  polytope.  Specification  of  the  coefficients  of  the  hydrodynamic  kets  then 
determine  the  cross  section  of  allowed  collision  outcomes,  or  the  kinetic  polytope,  parametrized  by 
the  coefficients  of  the  kinetic  kets  an.  Thus,  for  a  given  incoming  state,  h  can  be  thought  of  as  a 
scalar  function  on  the  kinetic  polytope.  We  denote  this  by  h(aT]),  and  we  demand  that  the  collision 
process  increase  this  function: 

Collision  Property  3  To  ensure  the  existence  and  stability  of  an  equilibrium  state,  a  collision  at 
site  x  must  not  decrease  the  restriction  of  the  function  h  to  the  kinetic  polytope. 

From  a  more  geometric  point  of  view,  note  that  the  (hyper)surfaces  of  constant  h(an)  provide 
a  codimension-one  foliation  of  the  kinetic  polytope.  These  (hyper)surfaces  degenerate  to  a  point 
where  h  reaches  a  maximum.  An  incoming  state  lies  on  such  a  (hyper)surface.  A  legitimate  collision 
is  required  to  map  such  an  incoming  state  to  an  outgoing  one  that  lies  inside,  or  at  least  on,  this 
(hyper)surface.  Clearly,  the  point  with  maximal  h  must  be  mapped  to  itself. 

Indeed,  we  can  subsume  Collision  Property  1  into  Collision  Property  3  by  constructing  the 
function  h  so  that  it  takes  a  minimum  value  on  the  boundary  of  the  polytope,  and  increases  to  a 
single  maximum  somewhere  inside.  The  master  polytope  is  defined  as  the  region  for  which  all  of 
the  Nj  (x)  ’s  are  nonnegative.  The  boundary  of  the  master  polytope  is  therefore  the  place  where  at 
least  one  of  the  Nj(x)'s  vanishes.  It  follows  that  Ylj  Nj(x)  is  constant  -  in  fact,  it  is  zero  -  on  the 
polytope  boundary.  More  generally,  if  (j(x)  for  j  =  1, . . .,  n  are  such  that 

0(0)  =  o,  (22) 

then  0  (Aj(x))  also  goes  to  zero  on  the  polytope  boundary.  To  ensure  that  there  is  only  one 
maximum  inside  the  convex  (master  or  kinetic)  polytope,  we  also  require  that  the  0  be  strictly 
increasing, 

Cj(x)  >  0,  (23) 

for  nonnegative  x.  Indeed,  from  Eqs.  (22)  and  (23)  it  follows  that 

Cj(x)  >  0,  (24) 

for  nonnegative  x.  Thus,  the  simplest  choice  would  be  to  make  h  (N(x))  a  function  of  []"  Cj  (Aj-(x)). 
To  make  this  consistent  with  Eq.  (21),  however,  we  see  that  the  functions  9j  should  be  chosen  to 
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be  the  logarithms  of  the  functions  (j,  so 


ft(N(x))  =  X:in[CiW(x))].  (25) 

3= 1 

This  goes  to  negative  infinity  on  the  boundaries  of  the  master  (and  hence  kinetic)  poly  topes,  and 
it  has  a  unique  maximum  in  the  interior. 

Collision  Property  4  A  valid  functional  form  for  h  is  given  by  Eq.  (25),  where  the  functions  (j 
obey  Eq.  (22),  Eq.  (23),  and  Eq.  (24).  . 

This  may  well  be  the  most  controversial  of  the  four  collision  properties  that  we  have  presented. 
Indeed,  it  is  not  necessary  to  assure  that  H  have  a  single  maximum  within  the  polytope.  We  shall 
discuss  an  H  function  that  violates  Collision  Property  4  in  Subsection  5.5. 

The  equilibrium  distribution  is  the  point  within  the  kinetic  polytope  where  h  has  its  maximum 
value.  If  we  denote  the  kinetic  parameters  by  aj,  then  we  can  find  this  point  by  demanding  that 
the  gradient  of  h  vanish, 


dh  CjWdNj  ="Cj(Nj) 

da,  Ci  (Nj)  da,  (Nj)  ‘ 


(26) 


If  we  denote  the  equilibrium  distribution  function  by  JVjq,  it  follows  that  the  covector  with  com¬ 
ponents  /(j  (iVjq)  lies  in  the  hydrodynamic  subspace, 


q  0VT) 
o  («T) 


~  1  j  • 


<7  =  1 


(27) 


The  nc  coefficients  Qa  must  be  chosen  so  that  the  values  of  the  conserved  quantities  are  as  desired, 


(Ar  |  Neq)  =  At  (28) 

for  t  —  1, . . . ,  nc.  Eqs.  (27)  and  (28)  can  be  regarded  as  n  +  nc  equations  for  the  n  +  nc  unknowns, 
Np  (for  j  =  1, . . .  ,  n)  and  Qa  (for  a  —  1, . . .,  nc).  In  general,  these  equations  are  nonlinear  and  it 
is  not  always  possible  to  write  the  iVjq  in  closed  form. 

The  equilibrium  distribution  is  something  over  which  the  model  builder  would  like  to  retain  as 
much  control  as  possible,  since  it  is  often  used  to  tailor  the  form  of  the  resultant  hydrodynamic  equa¬ 
tions.  For  example,  a  judicious  choice  of  the  equilibrium  distribution  function  is  required  to  obtain 
Galilean  invariant  hydrodynamic  equations  for  lattice  Boltzmann  models  of  fluids.  Customarily, 
lattice  Boltzmann  model  builders  have  simply  dictated  the  form  of  the  equilibrium  distribution 
function.  To  the  extent  that  this  is  necessary,  however,  it  is  more  in  keeping  with  the  philosophy 
espoused  in  this  paper  to  (try  to)  dictate  the  form  of  the  Lyapunov  functional,  rather  than  that  of 
the  equilibrium  distribution.  The  challenge  is  to  find  a  set  of  functions  (j,  subject  to  the  require¬ 
ments  of  Eqs.  (22),  (23)  and  (24),  such  that  Eq.  (27)  yields  the  desired  equilibrium  distribution. 
We  shall  return  to  this  problem  in  Subsection  5.5. 

While  this  formulation  of  an  entropic  principle  for  lattice  Boltzmann  models  seems  reasonable, 
note  that  it  is  rather  different  from  the  usual  one  encountered  in  kinetic  theory.  The  usual  choice 
there  would  be  that  of  Boltzmann,  h  =  £"=1  NjlnNj.  A  moment’s  examination  of  Eq.  (25) 
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indicates  that  this  would  correspond  to  the  choice  Cj(x)  =  %x but  this  function  decreases  for 
x  E  (0, 1/e)  and  hence  violates  Eq.  (23).  Thus,  while  it  is  possible  to  construct  lattice  Boltzmann 
models  with  Maxwell-Boltzmann  equilibria  [21],  the  Lyapunov  functionals  of  the  models  described 
in  this  paper  need  to  be  rather  different  from  those  commonly  used  in  kinetic  theory.  With  this 
caveat  in  mind,  we  shall  henceforth  abuse  terminology  by  taking  “H  function”  and  “Lyapunov 
functional”  to  be  synonymous. 

2.8  Collision  Operator 

Finally,  we  turn  our  attention  to  the  construction  of  a  collision  operator  that  is  in  keeping  with  the 
collision  properties  described  above.  Obviously,  there  are  many  ways  to  describe  maps  from  the 
kinetic  polytopes  into  (and  perhaps  onto)  themselves,  that  do  not  decrease  an  H  function.  For  the 
purposes  of  this  paper,  however,  we  restrict  our  attention  to  the  BGK  form  of  collision  operator  in 
which  the  outgoing  state  is  a  linear  combination  of  the  incoming  state  and  the  equilbrium, 

|C)  =  i(|Neq)-|N)),  (29) 

where  r  is  called  the  relaxation  time.  Usually  r  is  taken  to  be  constant,  but  more  generally  it  may 
depend  on  the  conserved  quantities,  and  most  generally  on  all  of  the  Nj. 

If  we  combine  Eqs.  (29)  and  (3),  we  get  the  BGK  equation, 

|N')  =  |N)  +  i(|Neq)-|N)).  (30) 

From  a  geometric  point  of  view,  this  equation  tells  us  to  draw  a  line  in  the  kinetic  polytope  from 
the  position  of  the  incoming  state  |N)  through  the  equilibrium  state  |Neq).  The  final  state  is  a 
weighted  combination  of  these  two  states  and  hence  lies  on  this  line.  The  incoming  state  is  weighted 
by  1  —  1/t,  and  the  equilibrium  state  is  weighted  by  1/r.  Thus,  for  r  >  1  the  outgoing  state  lies 
somewhere  on  the  segment  between  the  incoming  state  and  the  equilibrium.  Since  both  of  these 
states  are  in  the  kinetic  polytope,  and  since  this  polytope  is  convex,  the  outgoing  state  must  lie 
within  it  as  well.  Thus  Collision  Property  1  is  satisfied.  Moreover,  since  the  restriction  of  h  to  this 
segment  increases  as  one  approaches  the  equilibrium,  Collision  Property  3  is  also  satisfied. 

The  instabilities  associated  with  lattice  BGK  operators  arise  because  practitioners  try  to  over¬ 
relax  them.  A  large  class  of  lattice  BGK  models  for  the  Navier-Stokes  equations  have  shear  viscosity 
v  oc  (r  —  1/2).  In  an  effort  to  achieve  lower  viscosity  (and  hence  higher  Reynolds  number),  prac¬ 
titioners  set  t  to  values  between  1/2  and  unity.  In  this  situation,  the  outgoing  state  “overshoots” 
the  equilibrium  and  lies  on  the  other  side  of  the  polytope,  opposite  the  equilibrium  state.  For 
sufficiently  small  r,  the  outgoing  state  will  lie  on  a  contour  of  h  that  is  lower  than  that  of  the 
incoming  state,  thereby  violating  Collision  Property  3.  Still  smaller  values  of  r  might  cause  the 
outgoing  state  to  lie  outside  the  kinetic  poly  tope  altogether,  thereby  violating  Collision  Property  1. 
In  either  case,  numerical  instability  is  likely  to  result. 

2.9  Entropically  Stabilized  BGK  Operators 

A  potential  solution  to  this  problem  is  suggested  by  our  geometrical  viewpoint.  For  t  —  1  the 
outgoing  state  is  the  equilibrium,  for  which  h  is  maximal.  As  r  is  decreased  from  unity,  the  final 
value  of  h  decreases  from  its  maximal  value.  At  some  value  of  r  less  than  unity,  the  outgoing  value 
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of  h  will  be  equal  to  the  incoming  one.  In  order  to  respect  Collision  Property  3,  we  must  not  make 
r  any  lower  than  this  value,  given  by  the  solution  to  the  equation 

h  (N')  =  h  (N) ,  (31) 

where  |N')  is  given  by  Eq.  (30).  This  may  be  regarded  as  a  nonlinear  algebraic  equation  for  the 
single  scalar  unknown  r.  Indeed,  this  limitation  on  r  was  suggested  independently  by  Karlin  et 
al.  [22],  and  by  Chen  and  Teixeira  [23].  It  is  straightforward  to  find  the  solution  to  this  equation 
numerically,  since  it  is  easily  bounded:  We  know  that  the  desired  solution  has  an  upper  bound  of 
unity.  The  lower  bound  will  be  that  for  which  the  solution  leaves  the  kinetic  polytope;  this  happens 
when  one  of  the  N" s  vanish,  or  equivalently  when  r  is  equal  to  the  largest  value  of  1  —  Np/Nj, 
for  j  e  {1, . . n},  that  lies  between  zero  and  one.  Given  these  two  bounds  on  r,  the  regula  falsi 
algorithm  will  reliably  locate  the  desired  solution.  Call  this  solution  r*(N).  This  will  generally  be 
a  function  of  the  incoming  state.  A  useful  way  to  parametrize  r  is  then  to  write 

r(N)  =  ^29,  (32) 

where  0  <  k  <  1  is  a  constant  parameter.  The  case  k  — »  0  corresponds  to  r  — >  oo  so  that  the 
collision  operator  vanishes;  the  case  k  -*  1  corresponds  to  r  — >  r*,  which  is  the  largest  value 
possible  that  respects  Collision  Property  3.  Thus,  the  entropically  stabilized  version  of  the  lattice 
BGK  equation  is 

iN'>  =  |N>+^k)(|Neq)-|N))>  (33) 

where  r*  is  the  solution  to  Eq.  (31),  and  0  <  k  <  1. 

Of  course,  making  r  a  nontrivial  function  of  the  incoming  state  will  impact  the  hydrodynamic 
equations  derived  from  the  model.  The  challenge  to  the  designer  of  entropic  lattice  Boltzmann 
models  is  then  to  choose  the  Q  very  judiciously,  so  that  Eq.  (33)  yields  the  desired  hydrodynamic 
equations,  while  stability  is  guaranteed  by  keeping  0  <  k  <  1. 

In  constructing  a  lattice  Boltzmann  model  in  this  fashion,  as  opposed  to  using  the  simpler 
prescription  of  specifying  r,  it  may  be  argued  that  we  are  relinquishing  some  control  over  the 

transport  coefficients.  After  all,  if  v  oc  (r  —  1/2),  it  appears  that  we  can  specify  arbitrarily  sirihll 

viscosity  by  reducing  r.  In  fact,  this  is  not  the  case  since  uncontrolled  instabilities  are  known  to  set 
in  for  t  well  above  1  /2.  The  objective  of  entropic  lattice  Boltzmann  models  is  to  allow  the  user  some 
ability  to  overrelax  the  collision  process,  without  sacrificing  stability.  The  price  that  one  may  have 
to  pay  for  this  stable  overrelaxation  is  living  with  a  bounded  transport  coefficient.  Thus,  entropic 
lattice  Boltzmann  models  for  fluid  flow  may  be  restricted  to  minimum  values  of  viscosity.  As  we 
shall  show,  however,  for  certain  entropically  stabilized  lattice  Boltzmann  models,  this  minimum 
value  can  actually  be  zero.  This  allows  for  fully  explicit,  perfectly  conservative,  absolutely  stable 
algorithms  at  arbitrarily  small  transport  coefficient. 

Finally,  we  note  that  this  particular  prescription  is  only  one  way  of  creating  a  stable  lattice 
Boltzmann  algorithm.  In  fact,  any  map  obeying  Collision  Properties  1  and  3  will  work.  More 
general  mappings  of  polytopes  to  themselves  can  and  should  be  considered,  and  we  leave  this  to 
future  study. 
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3  One-Dimensional  Diffusion  Model 


Having  discussed  entropic  lattice  Boltzmann  models  in  general  terms,  we  now  apply  the  formalism  to 
mass  diffusion  in  one  dimension.  In  elementary  books  on  numerical  analysis,  it  is  demonstrated  that 
the  fully  explicit  finite-difference  approximation  to  the  one-dimensional  diffusion  equation  is  stable 
only  if  the  Courant  condition  [24],  At  <  (Ax)2/ {2 D),  is  satisfied;  note  that  this  places  an  upper 
bound  on  the  transport  coefficient.  It  is  also  shown  that  this  condition  may  be  removed  by  adopting 
an  implicit  differencing  scheme,  such  as  that  of  Crank  and  Nicolson  [25] ,  or  an  alternating-direction 
implicit  scheme,  such  as  that  of  Peaceman  and  Rachford  [26] .  The  DuFort-Frankel  algorithm  [27]  is 
fully  explicit  and  unconditionally  stable,  but  it  achieves  this  by  a  differencing  scheme  that  involves 
three  time  steps,  even  though  the  diffusion  equation  is  first-order  in  time. 

For  the  problem  of  achieving  high  Reynolds  number,  one  would  like  the  transport  coefficient  to 
be  as  small  as  possible.  We  shall  show  that  the  entropic  lattice  Boltzmann  algorithm  provides  a 
fully  explicit,  perfectly  conservative,  two-time-step  algorithm  that  is  absolutely  stable  for  arbitrarily 
small  transport  coefficient.  While  this  result  seems  significant  in  and  of  itself,  part  of  our  purpose  is 
pedagogical.  In  the  course  of  our  development  of  this  model,  we  shall  discuss  optimal  conservation 
representations  and  present  the  Fourier-Motzkin  method  for  visualizing  the  Piaster  and  kinetic 
polytopes.  Nongeneric  features  of  the  example  are  noted,  in  preparation  for  more  sophisticated 
examples  in  subsequent  sections. 


3.1  Description  of  Model 


We  suppose  that  we  have  a  regular  one-dimensional  lattice  C  whose  sites  x  6  C  are  occupied  by 
particles  whose  velocities  may  take  on  one  of  only  n  —  3  discrete  values,  namely  —1,0  and  +1.  We 
abuse  notation  slightly  by  letting  j  take  its  values  from  the  set  of  symbols  {—,0,+},  so  that  we 
may  write  the  single-particle  distribution  as  Nj(x).  Thus,  the  state  of  a  given  site  is  captured  by 
the  ket, 

(  N-  \ 


|N)  = 


N0 


(34) 


\ 


iV+ 


where  we  have  suppressed  the  dependence  on  the  coordinate  x  for  simplicity. 

Diffusion  conserves  mass,  so  we  suppose  that  the  mass  per  unit  lattice  site  (mass  density  in 
lattice  units), 

p  =  iV_  +  No  +  N+  =  (p  |  N) ,  (35) 


is  conserved  by  the  collisions.  Here  we  have  introduced  the  hydrodynamic  “bra,” 


(p|  -  (  +1  +1  +1  )  •  (36) 

(37) 

This  is  the  one  and  only  hydrodynamic  degree  of  freedom  in  this  example.  Because  there  are  a 
total  of  three  degrees  of  freedom,  the  other  two  must  be  kinetic  in  nature.  To  span  these  kinetic 
degrees  of  freedom  and  thereby  make  the  bra  basis  complete,  we  introduce  the  linearly  independent 
bras, 


(a|  =  (  -1  +2  -1  ) 

(38) 

01  =  (  +i  o  -l  ) . 

(39) 
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Next,  we  form  a  matrix  of  the  three  bras  and  invert  it  to  get  the  dual  basis  of  kets, 

-l  /  . .  .  -  .  \  -1 


0>  H  1/3)  )  = 


(  (P\ 

(a  j 

V  </3| 


From  Eq.  (40)  we  identify  the  linearly  independent  basis  kets 


l  p>  =  * 


+2  -1  +3 
+2  +2  0 
+2  -1  -3 


1/3)  =  5 


(40) 


(41) 


The  first  of  these  is  a  hydrodynamic  basis  ket,  while  the  last  two  are  kinetic  basis  kets. 

One  nongeneric  feature  should  be  noted:  In  this  example,  we  were  able  to  choose  (a|  and 
(/3|  so  that  each  ket  is  proportional  to  the  transpose  of  a  corresponding  bra.  Such  a  choice  is 
convenient  but  unnecessary.  The  only  real  requirement  in  choosing  the  kinetic  bras  is  that  they  be 
linearly  independent  of  each  other  and  of  the  hydrodynamic  bras.  The  next  example  will  illustrate 
a  situation  in  which  there  is  no  obvious  correspondence  between  the  individual  bras  and  kets. 

We  can  expand  the  state  ket  in  this  basis  as  follows 


|N>  =■  p\p)  +  a  |a)  +  /?  |/3) 


P 

3 


/ 

l 


i-!  + 

1  +  a 


(42) 


where  we  have  defined  a  =  a/ p  and  /3  =  P/p.  The  coefficients  p,  a  and  /?  (or  equivalently  p,  a  and 
/ 3 )  constitute  the  conservation  representation.  The  inverse  of  this  transformation  is  seen  to  be 


p  =  (p  |  N)  =  +AL  +  N0  +  N+ 

a  =  (a  |  N)  =  -N-  +  2iV0  -  N+  (43) 

/ 3  =  {/3  |  N)  =  +N-  -  N+. 


3.2  Nonnegativity 


The  nonnegativity  of  the  components  of  the  single-particle  distribution,  |N),  places  inequality 
constraints  on  the  parameters  p,  a  and  /3.  For  example,  it  is  clear  that  p  >  0.  Referring  to 
Eq.  (42),  we  see  that  we  must  also  demand 


0 

0 

0 


< 

"‘I 

< 

1  +  a 

< 

*-! 

+! 


2  ' 


(44) 


It  is  easy  to  see  that  all  three  of  these  inequalities  may  be  subsumed  by  the  single  statement 


-  1  <  a  <  2  -  3  \/3\ . 


(45) 


This  restricts  the  kinetic  parameters  to  a  triangular  region  in  the  (a,  (3)  plane,  as  is  shown  in  Fig.  1. 
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Polytope  for  D=1  Diffusion 


Figure  1:  The  Kinetic  Poly  tope:  Nonnegativity  of  the  distribution  function  requires  that  the 
kinetic  parameters  a  and  ft  lie  in  the  shaded  triangular  region. 


Note  that  the  effect  of  varying  p  at  constant  a  and  ft  is  to  simply  scale  the  components  of  the 
distribution  function.  The  triangular  region  bounding  the  parameters  a  and  ft  is  then  independent 
of  p.  That  is,  the  bounds  on  the  barred  kinetic  parameters  do  not  depend  on  the  hydrodynamic 
parameter,  so  there  is  really  no  need  for  distinction  between  the  master  and  kinetic  polytopes.  It 
should  be  noted  that  this  is  another  nongeneric  feature  of  the  present  model.  While  it  is  always 
possible  to  scale  out  by  a  single  nonnegative-definite  hydrodynamic  parameter,  as  we  have  done 
with  p  in  this  case,  more  sophisticated  models  will  have  several  hydrodynamic  parameters.  In 
such  cases,  the  shape  of  the  region  bounding  the  kinetic  parameters  will  depend  on  the  remaining 
hydrodynamic  parameters.  We  shall  see  this  more  clearly  in  the  example  of  Section  4. 

3.3  Optimality  of  Representation 

In  this  example,  there  is  clearly  a  great  deal  of  latitude  in  the  choice  of  the  kinetic  bras.  The  only 
requirement  that  we  have  placed  on  these  is  that  they  be  linearly  independent  of  the  hydrodynamic 
bra  and  of  each  other.  This  raises  the  question  of  whether  there  is  some  optimal  choice  that  might 
be  made  for  these.  Of  course,  this  depends  entirely  on  what  is  meant  by  “optimal”  in  this  context. 

It  has  been  noted  that  the  original  representation  of  |N)  is  more  natural  for  expressing  the 
constraint  of  nonnegativity,  while  the  conservation  representation  is  more  natural  for  expressing 
the  collision  process.  For  this  reason,  any  computer  implementation  of  entropic  lattice  Boltzmann 
methods  will  require  frequent  transformations  between  the  two  representations.  This  transforma¬ 
tion  is  precisely  what  we  have  worked  out  in  Eqs.  (42)  and  (43)  above.  Thus,  one  natural  figure  of 
merit  is  the  number  of  arithmetic  operations  required  to  perform  such  transformations. 

To  investigate  this  question,  rather  than  requiring  that  the  kinetic  bras  be  given  by  Eqs.  (38) 
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and  (39),  we  leave  them  in  the  general  form 


(<*|  =  ^  aq  c*2  as  J 

<01  =  (  01  02  03  )  • 

When  we  invert  to  get  the  kets,  the  results  are 


*  <*203  -  <*302  \ 

(  02  -  03  ^ 

1  1 

^  Q!3  —  0-2  N 

Ip)  =  i  | 

<*301  -  <*103  |<*)  =  S 

03-01 

H<1 

II 

ft 

a\  -  03 

^  <*102  -  <*201  j 

(  01-02  ) 

f  1 

^  a2  -  a\  y 

where  we  have  defined  the  determinant 


(46) 

(47) 


(48) 


A  =  <*1  (02  -  03)  +  <*2  (03  -  0i)  +  <*3  (/?!  -  02)  . 


(49) 


One  way  to  simplify  the  transformation  process  would  be  to  find  a  representation  in  which  as  many 
components  as  possible  of  the  kinetic  bras  and  kets  vanish,  without  making  the  transformation 
singular.  This  means  that  we  want  to  choose  the  aj's  and  fy's  such  that  A  ^  0,  while  making 
vanish  as  many  as  possible  of  the  following  twelve  quantities: 


<*1, 

01, 


as  -  <*2, 
02  -  03. 


<*2, 

02, 

<*i  -  ots , 


03  ~  01, 


<*3, 

03, 

oi2  ~  <*i , 


01  -  02- 


In  this  example,  it  is  straightforward  to  see  that  there  are  a  number  of  ways  to  make  six  of  these 
quantities  equal  to  zero.  For  example,  we  could  choose 


resulting  in  A  =  c*i/?2  7^  0  and 


0:3 

02 


<*1  7^  0, 

<*2  =  0, 

<*3  =  0, 

01  =  0, 

02  7^  0, 

03  =  0, 

(50) 

-<*2=0, 

<*1  -  <*3  7^  0, 

<*2  —  <*1  #  0, 

-  03  t  0, 

03  -  01  =  0, 

01  ~  02  ±  0- 

(51) 

This  choice  also  has  the  added  virtue  of  making  the  hydrodynamic  ket  equal  to 


(  0  \ 

0 

l1/ 


(52) 


which  also  has  two  vanishing  components. 

The  computation  of  the  |N)  components  from  the  parameters  p,  a  and  0  is  thus  reduced  to 


=  pa 

=  pfi 

=  p{\-a-fj), 


N- 

N0 

N+ 
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involving  a  total  of  two  multiplications  and  two  additions/subtractions;  this  may  be  contrasted 
with  Eqs.  (42)  which  involve  six  multiplications  and  five  additions/subtractions.  Likewise,  using 
the  choice  of  Eqs.  (50),  the  computation  of  the  parameters  p,  a  and  /3  from  the  N  is  reduced  to 

p  =  N_  +  N0  +  N+ 
a  =  N-/p 

P  =  N0/p, 


involving  two  additions/subtractions  and  two  divisions;  this  may  be  contrasted  to  Eqs.  (43)  which 
involve  five  additions/subtractions,  one  multiplication,  and  two  divisions. 

Clearly,  more  sophisticated  figures  of  merit  could  be  devised  to  optimize  the  transformations 
used  in  lattice  Boltzmann  computations.  As  we  have  noted,  vanishing  components  of  the  kinetic 
bras  and  kets  eliminate  the  addition/subtraction  of  terms.  Likewise,  components  equal  to  ±1  do 
require  an  addition/subtraction,  but  not  a  multiplication,  and  this  fact  could  be  taken  into  account 
in  a  more  refined  figure  of  merit.  The  computation  of  the  collision  outcome  is  the  principal  “inner 
loop”  of  a  lattice  Boltzmann  computation,  insofar  as  it  must  be  performed  at  each  site  of  a  spatial 
grid  at  each  time  step.  Such  considerations  may  be  especially  important  for  lattice  Boltzmann 
models  with  large  numbers  of  velocities. 


3.4  Fourier-Motzkin  Elimination 

For  this  simple  example,  we  had  no  difficulty  visualizing  the  master  and  kinetic  polytopes.  In 
preparation  for  the  succeeding  sections  where  the  task  will  be  substantially  more  difficult,  however, 
we  take  this  opportunity  to  introduce  the  Fourier-Motzkin  elimination  method  for  this  purpose. 
The  algorithm  consists  of  the  following  sequence  of  steps: 

1.  We  rewrite  the  set  of  inequalities,  Eq.  (44)  in  matrix  format, 


/-i 

+  3 
'2 

+1  \ 

( 

a  \ 

+i 

0 

+1 

P 

\ 

3 

2 

+1/ 

K 

w 

(53) 


We  have  adopted  the  convention  of  including  constant  terms  in  an  extra  column  of  the  matrix, 
using  the  device  of  appending  1  to  the  column  vector  of  unknowns.  In  general,  there  will  be 
m  inequalities  for  n  unknowns,  and  the  matrix  will  be  of  size  m  x  (n  +  1). 


2. 


We  scale  each  inequality  by  a  positive  factor  so  that  the  pivot  5  is  either  —1,  0  or  +1.  (Recall 
.  that  scaling  by  a  positive  factor  preserves  the  sense  of  an  inequality.)  We  then  reorder  the 
inequalities,  sorting  by  their  (scaled)  pivots  so  that  the  zero  pivots  are  last.  Beginning  with 
Eq.  (53),  this  yields 


/  -1  +3  +2  ^ 

(  «  \ 

-1  -3  +2 

P 

\  +1  0  +1  j 

K  1  ) 

(54) 


Note  that  there  were  no  zero  pivots  in  this  first  step.  We  did,  however,  reorder  the  inequalities 
so  that  those  with  pivot  —1  are  together,  and  preceed  that  with  pivot  +1. 

5As  in  discussions  of  Gaussian  elimination,  the  term  “pivot”  refers  to  the  first  nonzero  entry  in  a  row. 


21 


3.  We  now  add  all  pairs  of  inequalities,  such  that  the  first  member  of  the  pair  has  pivot  —1  and 
the  second  member  of  the  pair  has  pivot  +1 ,  and  we  append  the  new  zero-pivot  inequalities 
thus  obtained  to  the  system.  If  there  are  m_  inequalities  with  pivot  —  1,  and  m+  inequalities 
with  pivot  +1,  this  results  in  the  addition  of  m_m+  new  inequalities  to  the  system.  Since 
there  were  m  —  m_  —  m+  zero-pivot  inequalities  to  begin  with,  the  new  system  will  have  a 
total  of  mo  =  m  —  m_  —  m+  +  m_m+  zero-pivot  inequalities.  For  our  above  example,  m  —  3, 
m_  =  2  and  m+  =  1,  so  we  get  mo  =  2  zero-pivot  inequalities  in  the  system,  which  can  now 
be  written 

/  -1  +3  +2  \ 

-1-3+2  /  a  \ 

+1  0+1  (3  >0.  (55) 

0  +3  +3  ^  1  y 

\  o  -3  +3  / 

4.  Finally,  we  recurse  by  returning  to  step  2  for  the  mo  x  n  submatrix  obtained  by  tahing  only 
the  zero-pivot  inequalities  and  neglecting  the  first  unknown  (since  the  zero-pivot  inequalities 
do  not  involve  it  anyway).  We  continue  in  this  fashion  until  all  of  the  first  n  elements  of  the 
rows  of  the  matrix  are  zero. 

For  our  example,  the  recursion  step  asks  us  to  return  to  step  2  for  the  submatrix  indicated 
below 


As  before,  we  begin  by  normalizing  the  rows  and  sorting  them  by  pivot, 


In  this  example  there  is  now  one  row  of  the  submatrix  with  pivot  —1  and  one  with  pivot  +1,  so  we 
append  the  sum  of  these  to  the  system  to  get 


At  this  point,  the  first  n  elements  of  the  last  row  of  the  matrix  are  zero,  so  we  stop  the  recursion 
and  examine  the  last  row.  It  states  a  true  inequality,  namely  +2  >  0,  so  we  can  conclude  that  the 
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original  inequalities  are  not  mutually  exclusive;  that  is,  a  nonnull  polytope  of  solutions  will  exist. 
If  the  n  +  1  element  of  the  final  row  had  been  negative,  we  would  have  concluded  that  no  such 
solution  was  possible. 

Once  we  have  stopped  the  recursion  and  established  the  consistency  of  the  inequalities,  we 
can  get  the  final  set  of  inequalities  by  “back  substituting”  up  the  matrix.  The  last  two  nontrivial 
inequalities  tell  us  that  —  1  <  /3  <  +1,  or 

\P\  <  I-  (59) 

The  first  three  equations  then  yield  —  1  <  a  <  min  (2  —  3/3,  2  +  3/3)  which  is  equivalent  to  Eq.  (45), 
describing  the  triangular  region  of  Fig.  1. 

More  generally,  each  variable  will  have  a  lower  bound  that  is  the  maximum  of  all  the  lower 
bounds  determined  by  the  inequalities  with  pivot  +1,  and  an  upper  bound  that  is  the  minimum  of 
all  the  upper  bounds  determined  by  the  inequalities  with  pivot  —1.  The  arguments  of  the  maximum 
and  minimum  functions  will  depend  only  on  those  variables  that  have  already  been  bounded.  If 
we  were  using  the  technique  to  find  limits  of  integration  for  a  multiple  integral,  the  inequalities  at 
the  bottom  of  the  Fourier-Motzkin-eliminated  matrix  would  correspond  to  the  outermost  integrals. 
We  shall  revisit  this  technique  in  the  next  section. 


3.5  Collision  Operator 


To  illustrate  the  construction  of  an  entropically  stable  collision  operator  for  this  model,  we  adopt 
the  simplest  possible  H  function  by  taking  (j(x)  —  x.  Since  there  is  only  one  hydrodynamic  degree 
of  freedom,  Eqs.  (27)  reduce  to 

Ntq  =  /V0eq  =  A^q  -  i  (60) 

Q 

Eq.  (28)  then  tells  us  that  Q  —  3/ p,  so 


(61) 


Comparing  this  with  Eq.  (42),  we  see  that,  within  the  kinetic  polytope,  the  equilibrium  point  is 
the  origin,  aeq  =  /3eq  =  0.  A  contour  plot  of  the  H  function  on  the  kinetic  polytope  is  presented  in 
Fig.  2. 

The  lattice  BGK  equation  in  the  conservation  representation,  Eq.  (30),  is  then  equivalent  to 
the  prescription 

P'  =  P  (62) 

and 


where  2  =  1  —  1/r.  These  can  also  be  written  as  the  single  equation 
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Figure  2:  H  Function  on  the  Kinetic  Polytope:  Note  that  H  reaches  a  maximum  when 
a  =  /3  —  0,  and  goes  to  — oo  on  the  boundary  of  the  polytope. 


As  promised,  the  collision  is  dramatically  simplified  in  this  representation.  The  hydrodynamic  pa¬ 
rameter  is  unchanged;  only  the  two  kinetic  parameters  are  altered  by  the  collision.  Transformation 
back  to  the  original  representation  yields 
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or  equivalently 
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3.6  Entropic  Stabilization 

The  condition  for  marginal  entropic  stabilization,  Eq.  (31),  for  this  model  is  then 

In  N'+  +  In  Nq  +  In  N'_  =  lnN+  +  lnJV0  +  In  N_. 


(67) 


This  needs  to  be  solved  for  the  limiting  value  z  =  z+,  where  z*  —  1  —  1/r*.  Using  Eqs.  (42)  and 
(65),  this  becomes 


z*a  3  z*(3\  (  2*q  3  z*/3 

1  ~  ~2~  +  ~2~  )  I1  +  z*a)  !  1 


=  li-“  +  M 
1  2  2 


(1  +  d)  1- 
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3/5 


.  (68) 


This  appears  to  be  a  cubic  equation  for  z* ,  but  in  fact  z*  =  1  (corresponding  to  r»  — >  oo)  is  clearly 
a  root.  It  is  an  uninteresting  root  since  it  corresponds  to  no  collision.  Once  it  is  removed,  we  are 
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Figure  3:  Contours  of  z*(a, /?)  on  the  Kinetic  Polytope:  This  function  approaches  —1  when 
a  =  f3  =  0,  is  —  2  on  the  midpoints  of  the  sides  of  the  triangle,  and  is  —1/2  at  the  vertices. 


left  with  a  quadratic  for  z* ,  so  the  relevant  solution  is  easily  written  in  closed  form, 


3a2  —  a3  +  9/?2  +  9a/32  —  ^/3(a2  +  d3  +  3/3 2  —  9a/?2)  (3a2  —  a3  +  9/32  +  9  a/?2) 

________ 


(69) 


A  contour  plot  of  the  function  z*(d,  /?)  on  the  kinetic  polytope  is  displayed  in  Fig.  3.  There  it 
can  be  seen  that  z*  has  minimum  value  of  -2  at  the  midpoints  of  the  sides  of  the  triangle,  has 
maximum  value  of  —1/2  at  the  vertices,  and  approaches  —1  at  the  equilibrium  a  —  (3  =  0. 

The  entropically  stabilized  collision  operator  is  then  given  by  Eq.  (65),  with  r  — >  r*/«,  or 
equivalently 

z  — ►  1  —  k(1  —  z*)-  (70) 

It  is  useful  to  discuss  three  possible  choices  for  k: 


•  k  —>  0  means  that  z  — >  1  and  r  — >  oo.  This  is  the  uninteresting  limit  in  which  the  collision 
operator  vanishes  and  hydrodynamic  behavior  is  lost. 

•  k  =  1/2  means  that  z  =  (1  +  z*)/2  and  r  =  2/(1  —  z*).  This  is  interesting  because  near 
equilibrium,  where  z*  ~  —1,  it  coincides  with  the  prescription  r  =  1.  This  is  the  limit, 
mentioned  in  Subsection  2.9,  wherein  the  outgoing  state  is  the  local  equilibrium.  This  is  the 
smallest  value  of  r  for  which  the  standard  lattice  BGK  algorithm  is  guaranteed  to  be  stable. 
(In  fairness,  the  standard  algorithm  is  likely  to  be  stable  for  smaller  r;  it’s  just  that  this  is 
not  guaranteed.) 

•  k  — »  1  means  that  z  — »  z*  and  r  — »  r».  This  is  the  largest  value  of  r  for  which  the  entropic 
lattice  BGK  algorithm  is  guaranteed  to  be  stable. 
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3.7  Attainable  Transport  Coefficient 

The  complicated  form  of  Eq.  (69)  may  lead  one  to  believe  that  the  hydrodynamic  equation  (in 
this  case,  a  diffusion  equation)  would  be  very  difficult  to  derive  for  the  entropically  stabilized 
collision  operator.  In  fact,  this  is  not  the  case  at  all.  The  Chapman-Enskog  analysis  that  yields 
the  hydrodynamic  equations  requires  only  that  one  linearize  the  collision  operator  about  the  local 
equilibrium.  Since  the  local  equilibrium  has  a  =  j3  =  0,  it  is  clear  from  Eqs.  (63)  and  (64),  that 
only  the  value  of  2  at  equilibrium  will  enter.  Since  2*  — >  —1  at  the  equilibrium,  Eq.  (70)  indicates 
that  the  transport  coefficient  will  be  precisely  that  obtained  by  the  standard  lattice  BGK  algorithm 
with  r  =  l/(2/c). 

The  Chapman-Enskog  analysis  of  the  standard  lattice  BGK  algorithm  is  an  elementary  exer¬ 
cise  [28].  The  result  for  the  diffusivity  in  natural  lattice  units  is 


(71) 


Hence,  by  the  above  argument,  the  result  for  the  diffusivity  of  the  entropic  lattice  BGK  algorithm 
is 


(72) 


From  this,  we  can  clearly  see  the  benefit  of  the  entropic  algorithm.  The  standard  lattice  BGK 
algorithm  is  not  guaranteed  stable  for  r  <  1  or  k  >  1/2,  as  discussed  above.  From  Eq.  (71)  or 
(72),  respectively,  we  see  that  this  corresponds  to  a  minimum  diffusivity  of  1/6.  By  contrast,  the 
entropic  lattice  BGK  algorithm  loses  stability  only  when  «  =  1,  corresponding  to  D  —  0. 

It  is  remarkable  that  we  have  found  a  stable,  conservative,  explicit  algorithm  that  seems  to 
allow  the  transport  coefficient  to  become  arbitrarily  small.  After  all,  the  analysis  leading  to  the 
shear  viscosity  of  a  fluid  model  is  not  very  different  from  that  leading  to  the  diffusivity  above,  and 
arbitrarily  small  shear  viscosity  would  allow  for  arbitrarily  large  Reynolds  number.  In  the  field 
of  computational  fluid  dynamics,  however,  results  that  seem  too  good  to  be  true  usually  are  just 
that.  For  one  thing,  stability  does  not  imply  accuracy.  A  perfectly  stable  algorithm  is  not  terribly 
useful  if  it  does  not  converge  to  the  correct  answer.  Flows  at  higher  Reynolds  number  involve  ever 
smaller  eddies,  and  at  some  point  the  lattice  spacing  becomes  insufficient  to  resolve  these.  While 
it  is  comforting  to  have  an  algorithm  that  does  not  lose  stability  in  this  situation,  it  is  probably  a 
mistake  to  attach  much  physical  significance  to  its  results. 

Another  problem  is  that  the  time  required  for  the  system  to  come  to  equilibrium  goes  to  infinity 
as  k  — - >  1.  In  fact,  when  k  —  1  the  entropic  lattice  BGK  collision  operator  is  its  own  inverse:  If 
incoming  state  |N)  yields  outgoing  state  |N'),  then  incoming  state  |N')  would  yield  outgoing  state 
|N).  This  means  that  if  the  entire  lattice  were  initialized  slightly  away  from  equilibrium  with  no 
spatial  gradients  whatsoever,  it  would  simply  thrash  back  and  forth  between  two  states,  without 
ever  converging  to  the  desired  equilibrium.  When  the  time  required  to  achieve  local  equilibrium 
exceeds  time  scales  of  interest,  hydrodynamic  behavior  breaks  down. 

To  illustrate  this  problem,  we  simulated  the  entropically  stabilized  diffusion  model  with  the 
initial  density  profile 


p  =  Po  +  pi  sin 


(73) 


on  a  periodic  lattice.  We  used  lattice  size  L  =  32,  and  wavenumber  i  =  3  in  all  the  simulations 
reported  below.  We  initialized  the  state  of  each  site  in  a  local  equilibrium  N _  =  No  =  Ar+  =  p/3, 
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Figure  4:  Decay  of  sinusoidal  density  profile  for  k  =  0.9  (left)  and  k  —  0.99  (right). 
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Figure  5:  Decay  of  sinusoidal  density  profile  for  k  —  0.999 


without  the  Chapman-Enskog  correction  due  to  gradient.  Now  Eq.  (73)  is  a  solution  of  the  diffusion 
equation,  if  p\  decays  in  time  as  exp(— jt),  where 

So  we  fit  our  results  to  the  functional  form  of  Eq.  (73),  measured  pi(t),  and  did  a  least-squares 
fit  of  its  logarithm  to  a  linear  function  of  time  in  order  to  determine  7.  We  then  used  Eq.  (74)  to 
get  the  diffusion  coefficient,  and  we  compared  this  to  the  theoretical  value  provided  by  Eq.  (72)  for 
several  different  values  of  k,  approaching  unity. 

Fig.  4  shows  the  measured  value  of  pi(t)  for  k  =  0.9  and  k  =  0.99.  As  can  be  seen,  there  is 
an  initial  transient,  due  to  the  inadequacy  of  the  form  used  for  the  local  equilibrium  N-  =  No  = 
N+  =  p/3  in  the  presence  of  a  spatial  gradient.  The  left-hand  side  of  Fig.  5  shows  the  measured 
value  of  p\(t)  for  k  =  0.999,  and  the  right-hand  side  is  an  enlargement  of  the  transient  region. 
Fig.  6  shows  the  same  things  for  k  =  0.9999.  Note  that  the  transient  period  lengthens  considerably 
as  k  nears  unity.  We  could  have  substantially  reduced  these  transient  periods  had  we  used  the 
Chapman-Enskog  correction  to  the  local  equilibrium;  instead,  we  simply  waited  for  the  transient 
to  die  away  before  measuring  the  decay  constant  7. 

The  results  for  the  diffusivity  are  displayed  in  Table  1.  Note  that  we  find  reasonable  agreement 
until  we  get  to  the  cases  k  =  0.9999  and  k  —  0.99999.  To  see  what  is  going  wrong  for  these  cases, 
the  upper-left-hand  corner  of  Fig.  7  shows  the  measured  value  of  p\{t)  for  k  =  0.99999,  and  the 
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Figure  6:  Decay  of  sinusoidal  density  profile  for  k  —  0.9999 
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^theory 

■^meas 

0.9 

1.852  x  10~2 

1.75364  x  10“* 

0.99 

1.684  x  10~3 

1.67419  x  10“3 

0.999 

1.668  x  10~4 

1.66674  x  10“4 

0.9999 

1.667  x  10-5 

1.66198  x  10-5 

0.99999 

1.667  x  10"6 

2.28751  x  10~6 

Table  1:  Theoretical  and  measured  values  of  the  diffusivity  for  various  values  of  k. 


upper-right-hand  corner  shows  the  same  plot  with  a  reduced  range  for  the  ordinate.  In  addition 
to  the  oddly  shaped  envelope  of  the  initial  transient  period,  similar  to  that  seen  in  Fig.  6,  and 
enlarged  in  the  lower  part  of  this  figure,  we  note  that  the  plot  of  pi(t)  in  the  upper  right  seems 
to  be  increasing  in  thickness,  even  after  the  initial  transient.  Upon  further  enlargement,  shown  on 
the  left-hand  side  of  Fig.  8,  we  see  that  this  thickness  is  in  fact  a  subtle  transient  of  extremely  long 
duration.  The  right-hand  side  of  the  figure  shows  this  same  transient  at  the  longest  time  for  which 
the  simulation  in  Fig.  7  was  run.  Though  reduced  in  magnitude  (the  scales  of  the  ordinates  of  both 
graphs  in  Fig.  8  are  equal),  the  transient  is  still  present,  even  though  most  of  the  signal  itself  has 
diffused  away.  We  believe  that  this  transient  is  causing  the  anomaly  in  the  measured  diffusivity.  A 
much  longer  simulation  would  be  required  to  test  this  assertion. 

To  better  understand  this  loss  of  separation  between  the  kinetic  and  hydrodynamic  time  scales, 
note  that  the  time  required  for  the  signal  to  decay  away  is 


1  1 

^signal  CX  ^  OC  ^  OC 


K 

1  —  k' 


(75) 


The  time  required  for  the  transient  to  decay  to  1  /e  of  its  initial  value,  on  the  other  hand,  obeys 


^transient  =  —  ^ 

e 


(76) 


whence 


1 

^transient  i  • 

ln« 


(77) 


For  k  =  1  —  e,  where  e  is  small,  the  characteristic  times  in  Eqs.  (75)  and  (77)  both  go  like  1/e. 
Thus,  the  transients  begin  to  linger  for  hydrodynamic  time  scales,  and  it  becomes  impossible  to 
separate  the  kinetic  behavior  from  the  hydrodynamic  behavior. 
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When  the  time  required  for  the  decay  of  transients  is  comparable  to  hydrodynamic  time  scales 
of  interest,  the  usefulness  of  the  simulation  is  questionable.  Thus,  while  Table  1  indicates  good 
agreement  with  theory,  and  the  transport  coefficient  really  does  tend  to  zero  without  loss  of  stability, 
one  must  be  exceedingly  careful  in  the  preparation  of  the  initial  condition  to  exploit  this  feature  of 
entropic  lattice  Boltzmann  models.  In  particular,  transient  behavior  could  be  dramatically  reduced 
by  including  the  Chapman-Enskog  correction  to  first-order  in  the  gradient.  Indeed,  if  this  problem 
turns  out  to  be  the  only  obstacle  to  stable,  conservative,  explicit  algorithms  with  arbitrarily  small 
transport  coefficients,  higher-order  gradient  corrections  may  also  be  worth  investigating. 
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4  One-Dimensional  Compressible  Fluid  Model 


In  this  section,  we  apply  the  entropic  lattice  Boltzmann  method  to  a  simple  five-velocity  model  of 
fluid  dynamics  in  one  dimension,  first  considered  by  Renda  et  al.  in  1997  [14].  We  shall  find  that 
the  geometric  picture  is  much  richer  than  that  for  the  diffusion  model.  The  master  polytope  for 
this  model  is  four-dimensional,  and  we  shall  show  that  the  Fourier-Motzkin  algorithm  is  very  useful 
for  describing  it. 

4.1  Description  of  Model 

As  a  second  example,  we  consider  a  lattice  Boltzmann  model  for  a  one-dimensional  compressible 
fluid  with  conserved  mass,  momentum  and  energy,  first  studied  by  Renda  et  al.  [14].  The  velocity 
space  of  this  model  consists  of  five  discrete  values  of  velocity,  namely  0,  ±1,  and  ±2.  The  single¬ 
particle  distribution  at  site  x  with  velocity  j  E  {—2,  —1, 0,  +1,  +2}  is  denoted  by  Nj.  As  in  the  last 
example,  the  state  vector  (N)  at  a  given  site  x  will  be  denoted  by  a  ket, 


where  we  have  suppressed  the  dependence  on  x  for  simplicity. 

A  compressible  fluid  conserves  mass,  momentum  and  kinetic  energy,  so  we  suppose  that  the 
corresponding  densities 

P  =  E^  =  (P  IN)  (79) 

j=- 2 
+2 

*  =  E  i  N1  =  I  N>  (80) 

j=- 2 

e  =  E  4  =  <e  I  N>  ’  (81) 


are  conserved  by  the  collisions.  Here  we  have  introduced  the  bras 

H  =  i 

(*1  j  =  j 

(sly  =  i2/2, 


or  explicitly 


(p\  =  (  +1  +1  +1  +1  +1  ) 
<vr (  =  (  -2  -1  0  +1  +2  ) 
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(83) 

(84) 


(e|  =  ^  (  +4  +1  0  +1  +4  ) .  (85) 

These  are  the  three  hydrodynamic  degrees  of  freedom  in  this  example.  Because  there  are  a  total 
of  five  degrees  of  freedom,  the  other  two  must  be  kinetic  in  nature.  To  span  these  kinetic  degrees 
of  freedom  and  thereby  make  the  bra  basis  complete,  we  introduce  the  linearly  independent  bras, 

(«|  =  J  (  0  +1  0+1  0  )  (86) 

<01  =  K  0  -1  0+1  0).  (87) 

Next,  we  form  a  matrix  of  the  five  bras  and  invert  it  to  get  the  dual  basis  of  kets, 


(  I P)  k)  |e)  |a)  |/3)  )  = 


+2  +2  +2  +2  +2 
-4  -2  0  +2  +4 

+4  +1  0  +1  +4 

0+1  0+1  0 
0-1  0+1  0 


0  -1  +1  -1  +2 

0  0  0  +4  -4 

+4  0-2-6  0 

0  0  0  +4  +4 

0  +1  +1  -1  -2 


Prom  Eq.  (89)  we  identify  the  linearly  independent  basis  kets,  consisting  of  the  hydrodynamic  basis 


lff)  =  3 


k>  =  \  -2 

0 


and  the  kinetic  basis  kets 


l«)  =  3  -6 


10)  =  * 


Unlike  the  previous  example,  there  is  no  obvious  correspondence  between  the  individual  bra  and  ket 
basis  elements.  This  is  because  there  is  no  natural  notion  of  a  metric  in  this  space.  Nevertheless, 
as  noted  in  Section  2,  it  is  always  possible  to  construct  an  entire  ket  basis  from  an  entire  bra  basis, 
and  that  is  what  we  have  done  here. 

We  can  then  expand  the  state  ket  in  this  basis  as  follows 

/  |(e  —  #  —  a  +  2/3)  \ 
a  —  /3 

|N)  =  p  |p)  +  7T  |7r)  +  e  |e)  +  a  |a)  +  /?  j/3)  =  p  1  -  |  (e  +  3a)  ,  (92) 

a  +  /? 

v  |  (e  +  #  —  a  -  2/3)  / 

where  we  have  continued  to  use  an  overbar  to  denote  conserved  quantities  per  unit  mass;  for 
example,  n  =  ir/p  is  the  hydrodynamic  velocity.  As  with  the  last  example,  this  representation  of 
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the  distribution  function  |N)  has  the  virtue  of  making  the  conserved  quantities  manifest.  In  this 
case,  we  have  the  inverse  transformation 

p  =  (p  |  N)  =  +N- 2  +  N-i  +  No  +  N+\  +  N+  2 

7T  =  (7T  I  N)  =  -2AL2  -  N-i  +  N+1  +  2N+2 

e  =  (e|N>  =  +2N-2  +  N-1/2  +  N+1/2  +  2N+2  (93) 

a  =  (a|N)  =  +iV_1/2  +  iV+i/2 

/?  =  I  N)  =  -ALj/2  +  N+i/2. 

4.2  Nonnegativity 

We  now  consider  the  shape  of  the  master  and  kinetic  polytopes  for  the  compressible  fluid  model 
described  in  the  last  subsection.  Here  we  shall  find  a  much  richer  structure  than  the  simple  triangu¬ 
lar  region  that  we  found  for  the  diffusive  example  in  the  last  section  above.  Referring  to  Eq.  (92), 
we  see  that  that  nonnegativity  of  the  components  of  the  single-particle  distribution  function  is 
guaranteed  by  the  set  of  inequalities  p>  0  and 

0  <  §  —  7f  —  a +  20 

0  <  a-0 

0  <  2  —  e  —  3a  (94) 

0  <  a  +  0 
0  <  £  +  ff  —  a -20. 

These  define  the  master  polytope  in  the  four  dimensional  (if,  £,  a,  0)  space.  If  we  fix  the  hydro- 
dynamic  parameters,  w  and  £,  the  above  still  constitute  a  set  of  linear  inequalities  specifying  the 
corresponding  kinetic  polytope  in  the  two  dimensional  (<S,  0)  space.  Geometrically,  these  are  the 
intersections  of  the  master  polytope  with  the  (hyper)planes  of  fixed  hydrodynamic  parameters  n 
and  £.  Thus,  though  the  mass  density  p  still  scales  out  of  the  distribution  function,  the  shape  of  the 
kinetic  polytopes  in  the  (a,  0)  plane  will  depend  on  the  hydrodynamic  velocity  fr  and  the  kinetic 
energy  per  unit  mass  e. 

It  is  not  particularly  easy  to  visualize  the  shape  of  the  kinetic  polytope  in  the  (<5,  /?)  plane 
for  given  hydrodynamic  parameters  ff  and  e.  In  fact,  the  task  is  tedious  enough  that  we  have 
relegated  it  to  Appendix  A,  which  takes  the  reader  on  a  detailed  tour  of  the  four-dimensional 
master  polytope  and  its  kinetic  polytope  cross  sections.  To  summarize  the  conclusions  of  that 
Appendix,  the  projection  of  the  master  polytope  on  the  (ff,  i)  plane  is  illustrated  in  Fig.  9.  Values 
of  (ff,  S)  outside  the  shaded  regions  are  not  possible.  The  shape  of  the  kinetic  polytope  is  then  a 
triangle  when  (ff,  e)  is  in  the  most  lightly  shaded  regions  of  Fig.  9,  a  quadrilateral  when  (ff,  e)  is 
in  the  intermediately  shaded  regions  of  Fig.  9,  and  a  pentagon  when  (ff,  e)  is  in  the  most  darkly 
shaded  region  of  Fig.  9.  Illustrations  and  detailed  descriptions  of  the  kinetic  polytopes  for  particular 
values  of  (ff,  e)  in  each  of  the  seven  distinct  regions  with  ff  >  0  are  provided  in  the  seven  figures 
of  Appendix  A;  the  chart  below  Fig.  9  indicates  which  figure  corresponds  to  each  of  its  regions. 
Finally,  since  the  inequalities  are  invariant  under  (if,/?)  — >  (— ff,  —  /?),  the  kinetic  polytopes  for  ff  <  0 
are  obtained  from  those  for  ff  >  0  by  simply  reflecting  in  /?. 
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Figure  9:  Bounds  on  ff  and  e:  The  regions  with  the  lightest  shading  correspond  to  triangular 
regions  in  the  (<5,  (3)  plane,  those  with  intermediate  shading  correspond  to  quadrilateral  regions  in 
the  (a,  /?)  plane,  and  that  with  the  darkest  shading  corresponds  to  pentagonal  regions  in  the  (a,  (3) 
plane. 


4.3  Fourier-Motzkin  Analysis  for  the  ID  Compressible  Fluid  Model 

To  compute  the  master  polytope  for  the  compressible  fluid  model,  we  cast  the  inequalities  of 
Eqs.  (94)  into  matrix  format  as  follows, 
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(95) 


Here  we  have  put  the  hydrodynamic  parameters  n  and  e  after  the  kinetic  parameters  a  and  (3  in 
the  column  vector  of  unknowns,  for  reasons  that  will  become  apparent  in  a  moment.  Performing 
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the  Fourier-Motzkin  algorithm  on  this  set  of  inequalities,  we  finally  arrive  at  the  23  inequalities, 
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The  last  of  these  is  a  consistency  condition  which  is  seen  to  be  satisfied.  The  five  inequalities  prior 
to  that  yield 

max  ^0,  — i,  —1,  -4^  <  E  <  2, 
or 

0  <  e  <  2,  (97) 

which  explains  the  extent  of  Fig.  9  along  the  ordinate.  The  six  inequalities  prior  to  that  yield 


max 


^—2,  -2ej  <  7r  <  min  ^+2,  +|  4-  ^e,  +2e^  , 


or 


|7r|  <  min  ^2,  |  2e^  , 


(98) 


which  explains  the  extent  of  Fig.  9  along  the  abscissa.  Thus,  we  have  bounded  the  projection  of 
the  master  polytope  in  the  hydrodynamic  parameter  space. 

Given  hydrodynamic  parameters  in  the  allowed  region  thus  determined,  the  first  11  inequalities 
of  Eq.  (96)  then  specify  the  kinetic  polytopes.  (This  is  why  we  put  the  hydrodynamic  parameters 
after  the  kinetic  ones  in  the  column  vector  of  unknowns.)  Continuing  our  back  substitution,  the 
bounds  on  /3  are 


max 


~£  +  7r>  —  <  0  <  min  ^ 


.  „  _  -■£  _  _  e  +  7T 

<  P  <  mm  + - ,  +e  +  7r ,  T  — 

3  3 


)' 


(99) 
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Finally,  the  bounds  on  a  are 


max  (/3,  -(3)  <  5  <  min  (+2/3  —  n  +  e,  -2(3  +  #  +  £■).  (100) 

Eqs.  (99)  and  (100)  specify  the  shape  of  the  kinetic  polytope  in  the  (a,  (3)  plane,  given  the  hydrody¬ 
namic  parameters  it  and  e.  In  fact,  they  provide  a  more  succinct  way  of  determining  and  expressing 
the  bounds  of  both  the  master  polytope  and  the  kinetic  polytopes  than  the  analysis  presented  in 
Appendix  A. 
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5  Two  and  Three-Dimensional  Fluid  Dynamics 

In  this  section,  we  apply  our  methodology  to  much  richer  models,  namely  those  most  commonly 
used  for  two  and  three-dimensional  Navier-Stokes  hydrodynamics.  Thought  the  Fourier-Motzkin 
method  cannot  easily  be  applied  to  the  latter,  we  are  able  to  present  a  valid  collision  representation 
for  both  examples.  We  indicate  how  a  computer  simulation  for  such  models  might  be  implemented. 


5.1  FHP  Model 


The  diffusion  and  compressible  fluid  models  that  we  considered  above  had  three  and  five  particle 
velocities,  respectively.  Such  models  are  of  academic  and  pedagogical  interest,  but  not  particularly 
useful  as  simulational  tools.  The  smallest  model  that  has  been  used  for  serious  computational  fluid 
dynamics  calculations  in  two  dimensions  is  called  the  FHP  model.  This  is  a  six- velocity  model  on 
a  triangular  grid  which  has  been  shown  to  yield  isotropic,  incompressible  Navier-Stokes  flow  in  the 
appropriate  scaling  limit  [1,  2].  The  six  velocities  are  given  by 

Cj  =  x  cos  j  -f  y  sin  ,  (101) 


for  0  <  j  <  5.  The  model  conserves  mass  with  density 

5 


II 

£ 

(102) 

and  momentum  with  density 

5 

7T  =  Y,NjCj, 
j=0 

(103) 

so  the  hydrodynamic  bras  are 

(P\  = 

(  +1  +1  +1  +1 

+1 

+1) 

(104) 

(7TX|  = 

(  +1  +2  “5  -1 

1 

2 

+1 ) 

(105) 

(*Vl  = 

( 0  +4  +ir  0 

2 

-f ) 

(106) 

and  a  reasonable  choice  for  the  kinetic  bras  might  be 

(«!  = 

( +i  -l  +i  -l 

+1 

-0 

(107) 

<01  = 

(  _1  +2  +2  _1 

+2 

+3  ) 

(108) 

(tI  = 

(o  +1  0 

-#  )■ 

(109) 

(no) 

The  kets  can  then  be  found  as  with  the  preceeding  examples;  the  hydrodynamic  kets  are 
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/  +2  \ 
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+1 
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1  \  1 
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+i 

5  \Vx)  “  6 
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+i 

-1 

V  +i  / 

V  +1  / 
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I  *y) 
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o  \ 

+i 

+i 

0 

-1 
-1  / 


(111) 


and  the  kinetic  kets  are 


by) 


& 

6 


(  +1  \  /  -2  \ 

-1  +1 

-  W  =  l  —2 

+1  +1 

V  -i  /  \  +i  / 

The  general  representation  of  the  state  of  the  FHP  model  is  then 


( 


0\ 
+1 
-1 
0 

+1 

V  -i  ) 


IN) 


P  |p)  +  TTx  kx)  +  %  ky)  +  a  |a)  +  0 1/3)  +  7 17) 

(  1  +  27rx  +  a  —  20  \ 

1  +  7rx  +  \/3%  —  a  +  /?  +  %/37 

1  -  7rx  +  \/3%  +  a  +  0  -  y/3j 


1  —  2nx  —  a  —  2/3 


(112) 


(113) 


1  —  7rx  -  \/3i Ty  +  a  +  /3  +  y/3j 
\  1  +  TTX  -  V3ny  -  a  +  0  -  -v/37  J 

In  addition,  the  factors  of  -\/3  could  be  absorbed  into  ny  and  7  to  further  simplify  the  final  expression. 


5.2  FCHC  Model 

Early  attempts  to  generalize  the  FHP  model  to  obtain  a  three-dimensional  lattice  model  of  the 
Navier-Stokes  equations  were  thwarted  by  the  observation  that  no  regular  three-dimensional  lattice 
would  yield  the  required  isotropy  -  as  the  triangular  lattice  did  in  two  dimensions.  The  problem 
was  solved  by  noting  that  a  lattice  with  the  required  isotropy  existed  in  four  dimensions,  and  its 
three-dimensional  projection  could  be  used  to  construct  an  isotropic  model  [3],  The  four  dimen¬ 
sional  lattice  is  called  the  Face-Centered  Hypercubic  (FCHC)  lattice.  It  is  a  self-dual  lattice  with 
coordination  number  24.  It  can  be  defined  as  all  sites  ( i,j,k,l )  on  a  Cartesian  lattice  such  that 
i  +  j  +  k  + 1  is  even.  The  24  lattice  vectors  -  or,  equivalently,  the  24  lattice  sites  that  neighbor  the 
origin  (0, 0, 0, 0)  -  are  enumerated  in  Table  2. 

To  project  these  lattice  vectors  to  three-dimensional  space,  we  simply  ignore  the  ficticious  l 
coordinate.  When  we  do  this,  we  note  that  6  of  the  resulting  three-dimensional  lattice  vectors  will 
have  two  preimages  from  the  original  set  of  24  FCHC  lattice  vectors.  Rather  than  maintain  these 
as  separate  three-dimensional  vectors  (as  is  done  in  lattice-gas  studies),  it  is  possible  to  combine 
them  in  lattice  Boltzmann  studies,  simply  weighting  contributions  from  those  directions  by  a  factor 
of  two.  In  fact,  this  weighting  factor  can  be  thought  of  as  a  direction-dependent  particle  mass  rrij, 
so  that  the  three-dimensional  projection  of  the  FCHC  model  is  as  presented  in  Table  3. 

The  model  is  required  to  conserve  mass  and  three  components  of  momentum,  so  that  the  four 
hydrodynamic  bras  are  given  by 


H- 

=  rrij 

(114) 

{^x\j 

=  mjCjx 

(115) 

<7rvl i 

=  miCiy 

(116) 

(Kz\j 

=  miciz- 

(117) 

These  are  shown  in  the  first  four  rows  of  Table  4  in  Appendix  B,  followed  by  a  choice  of  twenty 
linearly  independent  kinetic  bras.  The  corresponding  kets  are  shown  in  Table  5  in  Appendix  B, 
and  these  can  be  used  to  construct  a  general  analytic  form  for  the  state  ket,  |N). 


38 


1  ciw 

cix 

C3y 

ciz 

1 

0 

+1 

+1 

0 

2 

0 

+1 

-1 

0 

3 

0 

-1 

+1 

0 

4 

0 

-1 

-1 

0 

5 

+1 

0 

0 

+1 

6 

-1 

0 

0 

+1 

7 

+1 

0 

0 

-1 

8 

-1 

0 

0 

-1 

9 

0 

+1 

0 

4-1 

10 

0 

+1 

0 

-1 

11 

0 

-1 

0 

+1 

j 

12 

0 

-1 

0 

-1 

13 

+1 

0 

+1 

0 

14 

-1 

0 

+1 

0 

15 

+1 

0 

-1 

0 

16 

-1 

0 

-1 

0 

17 

0 

0 

+1 

+1 

18 

0 

0 

+1 

-1 

19 

0 

0 

-1 

+1 

20 
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+1 
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24 

-1 

-1 

0 

0 

Table  2:  Lattice  vectors  of  the  FCHC  model. 
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0 

Table  3:  Lattice  vectors  and  masses  of  the  three-dimensional  projection  of  the  FCHC  model. 


5.3  Attainable  Transport  Coefficient 

As  with  the  one-dimensional  diffusion  model,  we  expect  that  the  entropic  lattice  Boltzmann  model 
will  allow  us  to  reduce  the  transport  coefficient  of  the  model.  For  fluids,  this  would  mean  smaller 
shear  viscosity  and  hence  higher  Reynolds  number.  To  see  if  this  is  possible,  we  consider  the  FHP 
model  in  the  incompressible  limit.  In  the  incompressible  limit,  the  Mach  number  scales  like  the 
Knudsen  number,  so  the  Chapman-Enskog  method  involves  only  the  equilibrium  distribution  at 
zero  momentum.  Taking  £j(x)  =  x,  it  is  straightforward  to  extremize  h  to  find  the  equilibrium 
distribution.  This  may  be  done  perturbatively  in  Mach  (equivalently,  Knudsen)  number,  and  the 
result  to  tenth  order  is 


eq  _ 


a 

peq 


2ttx  -  3ir^j  nfj  (l  +  2t ?x  +  2tt£)  +  O  (tt10)  (118) 

*1  ~  -  67rx7rv  +  7iy)  + 

2  (tt®  -  -  5 TT^TT^  +  7T®)  +  (t tI  +  nfy  2  (3tt*  -I-  -  5nfj  +  O  (tt10)  (119) 

2^1%  1  -  2  (ir*  -  +  4  (n*  -  n*')  -  2  ^  (3 tt*  +  j  +  O  (tt10)  .  (120) 


Thus  the  equilibrium  values  of  ft  and  7  are  of  second  order,  while  that  of  a  is  actually  third  order. 
For  the  computation  of  viscosity  in  the  incompressible  limit,  the  Chapman-Enskog  analysis  requires 
the  linearized  collision  operator  only  to  zeroth  order  in  the  Mach  number.  Thus,  for  the  purposes 
of  working  out  the  attainable  viscosity,  we  may  simply  set  the  equilibrium  values  of  all  three  kinetic 
parameters  to  zero.  It  follows  that  the  outgoing  kinetic  parameters  are  equal  to  the  incoming  ones 
multiplied  by  2  =  1  -  1/r. 
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We  can  now  construct  an  entropic  lattice  Boltzmann  model.  Using  Eq.  (113)  for  guidance,  we 
see  that  the  analog  of  Eq.  (68)  for  2*  is 

(l  +  2  Tiz  4 -a-  2/3)  (l  +  7fi  +  s/iHy  -  a  +  1 3+  v/37)  (l  -  nx  +  y/Sty  4  a  4  p  -  V37) 

(l  -  2nx  -  a-  2/3)  (l  -  Hx  -  ^/3nv  4  a  4  /3  4  v^y)  (l  4  nx  -  s/iHy  -  a  +  /3  -  = 

(l  4  2ftx  4  z*a  —  2 z*0)  (l  4  t*  4  V^TTy  —  z.a  4  z,/3  4  \/3z ,7)  (l  —  nx  4  s/3fiy  4  z»a  +  z,/3  — 

(l  —  2 ttx  —  z.a  —  2z,/3)  (l  —  7rx  —  \Z3ny  4  z.a  4  z.fi  4  %/3z,7)  (l  4  nx  -  \Z3tt y  —  z*a  4  z,j3  —  v/3z,7)  ,  (121) 

where  it  bears  repeating  that,  just  for  the  purposes  of  this  computation  of  viscosity,  we  have 
neglected  corrections  to  the  equilibrium  of  order  Mach  number  squared  6.  This  is  a  sixth-order 
equation  for  2*;  since  2,  =  1  is  clearly  a  root,  it  can  be  reduced  to  fifth  order.  In  spite  of  the  fact 
that  fifth  order  equations  generally  do  not  admit  solutions  in  radicals,  this  one  does  happen  to  do 
so.  The  result,  however,  is  complicated  and  uninspiring,  and  so  we  omit  it  here.  The  interested 
reader  can  simply  type  the  above  equation  into,  e.g.,  Mathematical© ,  and  use  the  Solve  utility  to 
see  the  result.  The  important  point  is  that  the  relevant  root  for  2,  approaches  —1  as  <5,/3, 7  — >  0. 

The  rest  of  the  analysis  is  exactly  as  it  was  for  the  one-dimensional  diffusion  model.  We  set 
2=1—  k(1  —  2*),  where  0  <  k  <  1.  Since  the  usual  lattice  Boltzmann  version  of  the  FHP  model 
has  the  shear  viscosity  proportional  to  r  —  1/2,  the  entropic  version  will  have  it  proportional  to 
1/k  -  1  which  goes  to  zero  as  k  approaches  unity.  Moreover,  absolute  stability  is  maintained  as 
this  limit  is  approached,  even  though  the  algorithm  is  fully  explicit  and  perfectly  conservative.  As 
with  the  diffusion  model,  we  are  prevented  from  achieving  arbitrarily  small  transport  coefficient 
only  by  considerations  of  accuracy  (resolution  of  the  smallest  eddies),  and  of  the  time  scale  for  the 
approach  to  equilibrium. 

5.4  Computer  Implementation 

We  pause  to  consider  the  computer  implementation  of  the  entropic  lattice  Boltzmann  algorithm 
for  the  FHP  fluid.  The  propagation  step  can  be  handled  in  precisely  the  same  fashion  that  it  is  for 
conventional  lattice  Boltzmann  simulations.  The  principal  difference  lies  in  the  implementation  of 
the  collision  operator,  and  that  is  what  we  describe  here.  Unlike  the  analysis  in  the  last  subsection, 
a  complete  computer  implementation  needs  to  solve  for  the  equilibrium  exactly  -  not  just  to  within 
the  Mach  number  squared. 

Given  the  six  quantities,  Nj  for  j  =  0,.  ..,5,  entering  a  site,  we  first  contract  them  with  the 
bras,  Eqs.  (104)  through  (109),  to  obtain  the  conservation  representation,  namely  p,  irx,  ny,  a,  /? 
and  7.  Given  ttx  and  ny,  it  is  then  necessary  to  solve  for  the  equilibrium  values  of  a,  ft  and  7. 
We  can  use  Eqs.  (118)  through  (120)  to  get  excellent  approximations  for  these,  and  then  use  a 
Newton-Raphson  iteration  to  get  the  correct  values  to  machine  precision.  Even  with  an  excellent 
initial  guess,  there  is  some  potential  for  loss  of  convergence  in  any  Newton-Raphson  iteration,  and 
so  it  might  be  worthwhile  to  investigate  better  ways  of  solving  for  these  values.  Assume  that  we 
can  do  this,  and  denote  the  equilibrium  values  with  an  “eq”  superscript. 

Given  the  equilibrium  values,  we  set 

a'  =  a  4  -  (aeq  -  a)  (122) 

_  T 

6Note  that  these  terms  of  order  Mach  number  squared  are  necessary  to  derive  the  remainder  of  the  Navier-Stokes 
equations. 
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=  ^  +  I(^q-^) 

T 

(123) 

i 

=  7  +  -  (7eq  -  7) , 

T 

(124) 

and  find  the  value  of  r  that  solves  the  nonlinear  equation  h(N')  =  h(N).  This  equation  will  look 
much  like  Eq.  (121),  except  that  <5eq,  /3eq  and  7eq  will  enter,  corresponding  to  corrections  of  order 
Mach  number  squared.  Since  this  nonlinear  equation  has  only  one  scalar  unknown,  namely  r,  it 
can  be  solved  using  an  iteration  method  that  is  guaranteed  to  converge,  such  as  regula  falsi.  We 
know  that  r  =  1  is  one  bound  on  the  solution;  the  other  bound  can  be  taken  as  the  point  at  which 
the  line  of  outgoing  states  intersects  the  boundary  of  the  kinetic  polytope.  This  means  finding  the 
largest  value  of  r  between  zero  and  one  for  which  Nj  =  0  for  some  j  e  {0, . . . ,  5).  Using  these  two 
bounds,  we  iterate  the  regula  falsi  algorithm  to  obtain  the  solution  t*.  We  then  set  r  =  t*/k,  and 
use  Eqs.  (122)  through  (124)  to  get  the  outgoing  values  of  the  kinetic  parameters.  We  finish  the 
collision  procedure  by  using  Eq.  (113)  to  get  the  outgoing  state  in  the  original  representation. 

5.5  Galilean  Invariance 

As  noted  in  the  Introduction,  practitioners  of  lattice  BGK  models  have  long  known  that  it  is 
necessary  to  tailor  the  equilibrium  distribution  function  in  a  certain  way  in  order  to  maintain 
Galilean  invariance.  Recall  that  the  zeroth  moment  (with  respect  to  the  velocity  vector)  of  the 
distribution  function  is  determined  by  the  mass  density.  Likewise,  the  first  moment  is  determined 
by  the  momentum  density.  For  a  lattice  BGK  model  of  an  ideal  gas  [29,  30],  the  trace  of  the  second 
moment  is  determined  by  the  kinetic  energy  density.  To  obtain  the  correct,  Galilean-invariant  form 
of  the  compressible  Navier-Stokes  equations,  it  is  also  necessary  [31]  to  mandate  the  rest  of  the 
second  moment,  the  third  moment,  and  the  trace  of  the  fourth  moment. 

While  it  is  possible  to  construct  equilibrium  distributions  with  these  moments  for  various  lat¬ 
tices,  no  lattice  model  has  ever  been  found  for  which  these  equilibria  are  guaranteed  to  be  nonneg¬ 
ative.  This  means  that  the  equilibrium  toward  which  we  would  like  to  relax  in  order  to  maintain 
Galilean  invariance  may  itself  lie  outside  of  the  master  polytope!  There  is  thus  no  way  that  such 
an  equilibrium  could  be  the  maximum  of  an  entropy  function  that  goes  to  — oo  on  the  polytope 
boundary,  such  as  do  the  ones  that  we  have  been  considering. 

In  such  a  situation,  we  may  ask  whether  it  is  possible  to  jettison  the  requirement  of  nonnegativ¬ 
ity,  while  still  demanding  absolute  stability.  In  fact,  our  formalism  suggests  an  interesting  way  of 
doing  just  that.  If  we  no  longer  demand  that  the  function  h  go  to  negative  infinity  on  the  polytope 
boundary,  then  we  need  no  longer  demand  Eq.  (22),  though  Eq.  (23)  is  still  useful  to  ensure  that 
the  equilibrium  is  unique,  and  Eq.  (24)  is  useful  to  keep  h  real- valued  in  light  of  Eq.  (25).  Thus, 
we  are  free  to  make  the  choice  Q(x)  is  exp[:r2/(2<?j)],  where  the  g/s  are  positive  constants.  Indeed, 
this  choice  is  similar  to  one  made  by  Karlin  et  al.  [32],  Eq.  (27)  is  then  the  set  of  n  equations, 

=  |\r>,,  (125) 

<7=1 

and  Eqs.  (28)  give  nc  additional  requirements.  One  may  then  try  to  solve  these  n  +  nc  equations 
for  the  n  +  nc  unknowns,  Q„  and  g-j,  subject  to  the  requirement  that  gj  >  0.  In  fact,  not  all  of 
these  unknowns  are  independent;  the  scaling  gj  — >  gjjx  and  Qa  — >  Qax,  where  x  is  arbitrary,  leaves 
Eq.  (125)  invariant,  so  there  are  more  requirements  than  equations,  and  hence  not  all  equilibria  are 
derivable  from  an  h  function  of  this  form.  Nevertheless,  many  interesting  equilibria  are  derivable 
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in  this  way,  and  this  approach  leads  immediately  to  the  construction  of  an  entropically  stabilized 
algorithm  for  such  equilibria.  Investigation  of  this  approach  is  work  in  progress  [33].  Its  successful 
resolution  may  yield  an  absolutely  stable  entropic  lattice  Boltzmann  algorithm  for  compressible 
flow,  albeit  one  that  allows  for  negative  Nj. 
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6  Conclusions 


In  summary,  we  have  presented  a  general  methodology  for  constructing  lattice  Boltzmann  models 
of  fluid  dynamics  that  conserve  a  given  set  of  moments  of  the  distribution  function.  The  method 
guarantees  the  nonnegativity  of  the  distribution  function,  and  allows  for  convenient  visualization  of 
the  dynamics  by  construction  of  the  polytope  of  allowed  states  using  Fourier-Motzkin  elimination. 
We  have  shown  how  such  models  can  be  endowed  with  a  Lyapunov  functional,  resulting  in  uncondi¬ 
tional  numerical  stability,  and  we  presented  various  choices  for  the  H  function  for  this  purpose.  We 
also  described  the  computer  implementation  of  such  models  in  some  detail.  For  lattice  Boltzmann 
models  of  diffusion  and  of  incompressible  fluid  flow,  we  were  able  to  present  fully  explicit,  perfectly 
conservative,  absolutely  stable  algorithms  that  appear  to  work  for  arbitrarily  small  transport  coef¬ 
ficient.  Indeed,  we  showed  that  the  limitations  on  attainable  transport  coefficients  for  these  models 
arise  from  considerations  of  accuracy,  rather  than  stability. 

There  are  numerous  avenues  for  extension  of  the  current  work.  We  close  by  presenting  a  list  of 
these.  We  hope  that  the  richness  and  promise  of  this  style  of  model  will  inspire  others  to  take  up 
some  of  these  questions. 

(i)  We  have  not  discussed  the  question  of  boundary  conditions  in  this  paper.  For  a  viscous  fluid, 
for  example,  collisions  at  a  solid  boundary  wall  must  conserve  mass,  but  not  momentum  or 
energy.  The  latter  two  quantities  can  enter  or  leave  the  domain  through  the  wall,  resulting 
in  drag,  lift,  surface  heating,  etc.  This  results  in  very  different  polytope  structure  at  the 
boundaries  than  in  the  interior.  In  addition,  there  is  nothing  preventing  entropy  from  entering 
or  leaving  through  the  wall,  so  the  entire  question  of  how  to  construct  h  functions  at  the 
boundary  sites  needs  to  be  revisited. 

(ii)  While  we  touched  on  the  question  of  optimality  of  conservation  representations,  clearly  much 
more  work  could  be  done  along  these  lines.  In  addition  to  making  it  computationally  easier 
to  transform  back  and  forth  between  the  original  and  conservation  representations,  one  might 
like  to  try  to  optimize  the  form  of  the  entropically  stabilized  collision  operator. 

(iii)  We  mentioned  reversible  models,  but  did  not  discuss  or  study  them  in  any  detail.  Our 
entropically  stabilized  models  are  reversible  in  the  limit  as  «  — ►  1,  but  this  is  precisely  the 
limit  in  which  they  begin  to  fail  to  simulate  the  parabolic  equations  of  interest,  namely  the 
diffusion  and  Navier-Stokes  equations.  It  would  be  interesting  to  know  if  such  reversible 
models  are  capable  of  simuating  time-symmetric  partial  differential  equations,  such  as  Euler’s 
equations  of  inviscid  fluid  dynamics. 

(iv)  The  nature  of  the  anomalies  encountered  as  k  — >  1  remains  to  be  completely  elucidated.  For 
example,  the  lower  right  plot  in  Fig.  7  appears  to  indicate  that  the  approach  to  equilibrium  is 
extremely  slow  for  this  model.  Does  it  eventually  reach  an  exponentially  decaying  equilibrium 
state  for  which  the  diffusion  coefficient  matches  that  of  the  theory?  Does  roundoff  error  play  a 
role  in  this  regime?  These  open  questions  will  require  substantially  more  numerical  simulation 
than  that  presented  here. 

(v)  Given  that  the  anomalies  as  k  — >  1  are  due  to  the  slowness  of  the  approach  to  the  equilibrium 
state,  much  effort  could  be  expended  in  trying  to  start  the  simulation  in  as  close  to  a  global 
equilibrium  state  as  possible.  This  means  taking  spatial  gradient  corrections  into  account.  At 
a  minimum,  the  Chapman-Enskog  corrections  -  accurate  to  first-order  in  the  spatial  gradient 
-  could  be  incorporated  with  little  effort. 
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(vi)  The  question  of  how  to  deal  with  interaction  potentials  in  this  formalism  could  be  addressed. 
At  the  level  of  mean-field  theory,  this  would  require  the  introduction  of  conserved  quantities 
that  are  quadratic  in  the  distribution  function.  The  allowed  regions  would  then  no  longer 
be  polytopes,  convex,  or  even  simply  connected.  Moreover,  multiple  local  equilibria  within 
the  allowed  region  may  be  expected.  For  near-equilibrium  dynamics,  this  might  give  rise  to  a 
Ginzburg-Landau  model  with  multiple  free-energy  minima.  This  is  thus  a  difficult  extension 
of  the  current  theory,  but  one  for  which  the  rewards  would  be  great. 

(vii)  One  could  try  to  construct  an  adaptive  mesh  refinement  (AMR)  version  of  this  algorithm 
for  which  the  nonnegativity  of  the  physical  densities  is  the  refinement  criterion.  Using  such 
an  approach,  it  might  be  possible  to  relax  the  criterion  of  nonnegativity  of  the  distribution 
function.  That  is,  we  could  carry  negative  distribution  function  values,  but  insist  that  all 
physical  densities  be  positive.  Should  the  propagation  step  threaten  to  create  a  negative 
density,  we  would  try  to  refine  the  lattice  locally  until  that  is  no  longer  the  case. 

(viii)  Finally,  it  would  be  interesting  to  consider  the  potential  utility  of  this  algorithm  for  the  con¬ 
struction  of  “eddy-viscosity”  subgrid  models  of  turbulence.  As  noted  at  several  points  in  this 
paper,  while  entropic  lattice  Boltzmann  algorithms  may  be  stable  for  arbitrarily  small  trans¬ 
port  coefficient,  they  may  begin  to  lose  accuracy.  The  smallest  eddy  sizes  in  three  dimensional 
Navier-Stokes  turbulence,  for  example,  scale  as  Re-3/4,  where  the  Reynolds  number  Re  goes 
inversely  with  viscosity.  Thus,  at  sufficiently  small  viscosity,  the  eddies  will  be  smaller  than  any 
fixed  grid  spacing.  If  one’s  goal  is  an  ab  initio  simulation  of  the  Navier-Stokes  equations,  this 
circumstance  is  grounds  for  rejecting  the  entropic  lattice  Boltzmann  solution  as  unphysical. 
On  the  other  hand,  stable,  coarse-grained,  turbulent  flow  must  obey  the  same  conservation 
laws  as  the  underlying  Navier-Stokes  equations.  Since  the  entropic  lattice  Boltzmann  colli¬ 
sion  operators  are  the  most  general  that  respect  these  conservation  laws,  it  is  interesting  to 
speculate  as  to  whether  some  physical  significance  might  yet  be  attached  to  their  solutions, 
even  when  the  smallest  eddies  are  not  resolved  by  the  lattice. 


Acknowlegements 

This  work  was  presented  by  BMB  at  the  Seventh  International  Conference  on  the  Discrete  Simu¬ 
lation  of  Fluids  in  Oxford,  England  in  July,  1998.  During  the  late  stages  of  writing  this  paper,  we 
became  aware  that  two  other  groups  [22,  23]  were  independently  working  on  requirements  of  the 
form  of  Eq.  (31).  As  the  works  were  conducted  independently,  there  are  numerous  difference  in 
style,  formulation  and  emphasis,  and  we  believe  that  the  resulting  papers,  including  this  one,  are 
nicely  complementary  to  one  another. 

It  is  a  pleasure  to  acknowledge  helpful  conversations  with  Frank  Alexander,  Bastien  Chopard, 
Nicos  Martys,  Sauro  Succi  and  Julia  Yeomans.  BMB  was  supported  in  part  by  an  IPA  assignment 
agreement  with  the  Air  Force  Research  Laboratory,  Hanscom  AFB,  Massachusetts,  and  in  part 
by  the  United  States  Air  Force  Office  of  Scientific  Research  under  grant  number  F49620-95-1- 
0285.  BMB  and  PVC  would  like  to  thank  Schlumberger  Research  and  the  EPSRC’s  Collaborative 
Computing  Project  No.  5  for  making  possible  BMB’s  visit  to  Schlumberger  Cambridge  Research 
and  the  University  of  Oxford,  where  this  work  was  begun  in  early  1997.  PVC  is  grateful  to  Wolfson 
College  and  the  Department  of  Theoretical  Physics,  University  of  Oxford,  for  a  Visiting  Fellowship 
(1996  -  1999). 


45 


References 

[1]  U.  Frisch,  B.  Hasslacher,  and  Y.  Pomeau,  Phys.  Rev.  Lett.  56  1505  (1986). 

[2]  S.  Wolfram,  J.  Stat.  Phys.  45,  471  (1986). 

[3]  U.  Frisch  et  ah,  Complex  Syst.  1,  648  (1987). 

[4]  D.H.  Rothman  and  S.  Zaleski,  Lattice-Gas  Automata:  Simple  Models  of  Complex  Hydrody¬ 
namics. ,  (Cambridge  University  Press,  1997). 

[5]  Henon,  M.,  Complex  Systems  1  (1987)  763. 

[6]  See,  e.g.,  B.M.  Boghosian,  P.V.  Coveney,  and  A.N.  Emerton,  Proc.  R.  Soc.  Lond.  A,  452, 
1221  (1996);  and  B.M.  Boghosian,  P.V.  Coveney  and  P.J.  Love,  “Three-dimensional  lattice- 
gas  model  for  amphiphilic  fluid  dynamics”,  Proc.  R.  Soc.  London  A ,  in  press  (1999). 

[7]  F.  Higuera  and  J.  Jimenez,  Europhys.  Lett.  9  (1989)  663. 

[8]  F.  Higuera,  S.  Sued  and  R.  Benzi,,  Europhys.  Lett.  9  (1989)  345. 

[9]  G.  McNamara  and  G.  Zanetti,  Phys.  Rev.  Lett.  61  (1989)  2332. 

[10]  S.  Succi,  R.  Benzi,  F.  Higuera,  Physica  D  47  (1991)  219-230. 

[11]  R.  Benzi,  S.  Succi,  and  M.  Vergassola  Phys.  Rep.,  222  145  (1992). 

[12]  E.  Gross,  D.  Bhatnager,  M.  Krook,  Phys.  Rev.  94  (1954)  511. 

[13]  Y.-H.  Qian,  D.  d’Humieres,  P.  Lallemand,  Europhys.  Lett.  17  (1992)  479-484. 

[14]  A.  Renda,  G.  Bella,  S.  Succi,  I.V.  Karlin,  Europhys.  Lett.  41  (1998)  279-283. 

[15]  B.M.  Boghosian,  P.V.  Coveney,  Comp.  Phys.  Comm,  (to  appear,  2000). 

[16]  Shan,  X.,  Chen,  H.,  Phys.  Rev.  E  47  (1993)  1815-1819. 

[17]  Swift,  M.R.,  Osborn,  W.R.,  Yeomans,  J.M.,  Phys.  Rev.  Lett.  75  (1995)  830-833;  Swift,  M.R., 
Orlandini,  S.E.,  Osborn,  W.R.,  Yeomans,  J.M.,  Phys.  Rev.  E  54  (1996)  5041-5052. 

[18]  N.  Martys,  Int.  J.  Mod.  Phys.  C 10  (1999)  1367-1382. 

[19]  D.  d’Humieres,  AIAA  Rarefied  Gas  Dynamics:  Theory  and  Applications,  in  Progress  in  As¬ 
tronautics  and  Aeronautics  159  (1992)  450-458. 

[20]  M.  Schechter,  Amer.  Math.  Monthly  105  (March,  1998)  246-251. 

[21]  A.J.  Wagner,  Europhys.  Lett.  44  (1998)  144-149. 

[22]  I.V.  Karlin,  A.  Ferrante  and  H.C.  Ottinger,  Europhys.  Lett.  47  (1999)  182-188. 

[23]  H.  Chen  and  C.  Teixeira,  “H-theorem  and  Origins  of  Instability  in  Thermal  Lattice  Boltzmann 
Models,”  preprint  (1999). 

[24]  R.  Courant,  K.  Friedrichs,  H.  Lewy,  Math.  Ann.  100  (1928)  32. 


46 


[25]  J.  Crank,  P.  Nicolson,  Proc.  Camb.  Phil.  Soc.  43  (1947)  50-67. 

[26]  D.W.  Peaceman,  H.H.  Rachford  Jr.,  J.  Soc.  Ind.  Appl.  Math.  3  (1955)  28. 

[27]  E.C.  DuFort,  S.P.  Frankel,  Math.  Tables  Aids  Comput.  7  (1953)  135-152. 

[28]  See,  for  example,  B.M.  Boghosian,  “The  Chapman-Enskog  Method  for  Lattice  Gases,”  in 
1989  Lectures  in  Complex  Systems,  Santa  Fe  Institute,  E.  Jen,  ed.  (1989). 

[29]  G.R.  McNamara,  A.L.  Garcia,  B.J.  Alder,  J.  Stat.  Phys.  87  (1997)  1111-1121. 

[30]  F.J.  Alexander,  A.L.  Garcia,  B.J.  Alder,  Phys.  Rev.  Lett.  74  (1995)  5212-5215. 

[31]  B.M.  Boghosian,  P.V.  Coveney,  Int.  J.  Mod.  Phys.  9  (1998)  1231-1246. 

[32]  I.  Karlin,  Phys.  Rev.  Lett.  81  (1998)  6. 

[33]  B.M.  Boghosian,  S.  Succi,  work  in  progress. 

1.1,  5.1  1.1,  5.1  1.1,  5.2  1.1  1  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.5,  4,  4.1  3  3  3  3  4  2.4  2.7,  2.7 
2.9,  6  2.9,  6  3  3  3  3  3.7  5.5  5.5  5.5  5.5  5.5 


47 


A  A  Tour  of  the  Master  Polytope  for  the  Compressible-Fluid 
Model 

The  five  inequalities  of  Eq.  (94)  can  be  recast  as 


0  <  |/3|  <  a  <  min 


( 


2  —  e 
3 


(126) 


This  implies 


e  >  a  +  1 2/3  —  7r|  >  a  >  |/3|  >  0 


(127) 


and 

£<2-35  <2  —  3|/3]  <  2,  (128) 

and  hence  a  >  0  and  £  €  [0, 2].  Moreover,  we  note  that  all  of  the  inequalities  are  invariant  under 
the  transformation 


7 r  — >  — 7T 

p  -  -p, 

so  that  we  may  restrict  our  attention  to  ff  >  0  without  loss  of  generality. 

Next  note  that  the  second  argument  of  the  min  function  in  Eq.  (126)  will  be  smaller  than  the 
first  if  0  <  e  <  (2  —  i)/ 3,  or  0  <  e  <  1/2.  In  this  case,  for  sufficiently  small  ff,  the  boundary  of  the 
polytope  in  the  (a,  /3)  plane  will  be  a  quadrilateral,  such  as  that  in  Fig.  10,  whose  lower  bound  in 
&  (shown  in  black)  is  |/3|  and  whose  upper  bound  in  a  (shown  in  gray)  is  £  —  |2/3  -  ff|.  The  allowed 
region  of  the  (a,  ft)  plane  is  the  shaded  area  in  between  these  bounds.  The  bottom  vertex  7  (a)  of 
this  quadrilateral  is  at  the  origin,  (a,/3)  =  (0,0),  and  the  upper  vertex  (c)  is  at  (a,j3)  =  (£,  ff/2). 
The  right  vertex  ( b )  is  then  at  /3  =  e  —  (2/3  — ff)  or  (a,/3)  =  ((£-fff)/3,  (£+ff)/3),  and  the  left  vertex 
( d )  is  then  at  —0  =  £  +  (2/3  -  f)  or  (a,  0)  =  ((e  -  7f ) /3,  (-e  +  7r)/3).  These  results  are  summarized 
in  the  table  included  in  Fig.  10. 

We  now  ask  for  what  range  of  £  and  n  the  above-described  quadrilateral  boundary  in  the  (a,  /3) 
plane  is  valid.  We  note  that  vertices  (a)  and  (6)  will  be  degenerate  when  e  =  —  7f,  and  that  vertices 
(a)  and  ( d )  will  be  degenerate  when  £  —  +7t.  Also,  it  is  easy  to  see  that  one  or  the  other  of 
these  two  degeneracies  will  occur  before  any  degeneracy  involving  vertex  (c).  It  follows  that  the 
above-described  quadrilateral  boundary  in  the  (a,/3)  plane  is  valid  only  for  |7f|  <  e  <  1/2.  This 
region  of  the  (ft,  e)  plane  is  the  shaded  triangle  AIH  in  Fig.  9. 

Beginning  in  triangle  AIH ,  boundary  AI  is  encountered  when  e  =  ft  so  that  vertices  (a)  and 
(d)  are  degenerate.  If  we  cross  this  boundary  the  allowed  region  in  the  (a,  /3)  plane  is  no  longer 
a  quadrilateral,  but  rather  becomes  a  triangle,  such  as  that  in  Fig.  11,  whose  upper  bound  in  a 
(shown  in  gray)  is  e  —  |2/3  —  ft|  and  whose  lower  bound  in  a  (shown  in  black)  is  /3.  The  bottom 
vertex  (e)  of  this  triangle  is  at  a  =  /3  =  i  +  (2/3  —  ft)  or  (a,  /3)  —  (ft  —  e,  ft  —  e);  this  vertex  replaces 
the  degenerate  vertices  (a)  and  ( d )  of  the  above-described  quadrilateral.  The  expressions  for  the 
coordinates  of  the  triangle’s  upper  vertex  (c)  and  its  right  vertex  ( b )  are  identical  to  those  of  the 
corresponding  vertices  of  the  quadrilateral.  These  results  are  summarized  in  the  table  included  in 
Fig.  11. 

Again,  we  ask  for  what  range  of  £  and  fr  the  above-described  triangular  boundary  in  the  (a,  (3) 
plane  is  valid.  We  note  that  vertices  (6)  and  (c)  will  be  degenerate  when  (£  +  fr) /3  =  s  or  e  =  fr/2; 

7Throughout  this  appendix,  we  denote  vertices  in  parentheses,  so  as  not  to  confuse  them  with  other  variables. 
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for  smaller  values  of  E  the  set  of  allowed  points  in  the  (a,  0)  plane  is  empty.  We  also  note  that,  just 
as  for  the  quadrilateral  of  Fig.  10,  the  first  argument  of  the  min  function  in  Eq.  (126)  will  become 
the  determining  factor  when  E  >  1/2,  also  invalidating  the  above  argument.  Thus,  the  allowed 
region  in  the  (a,  0)  plane  will  be  a  triangle  with  vertices  described  in  Fig.  11  only  if  (E,  tt)  is  in  the 
shaded  triangle  ABI  in  Fig.  9.  Since  the  region  AEH  is  the  image  of  ABI  under  the  map  n  — >  — ff, 
it  follows  that  the  allowed  region  of  the  (a,  0)  plane  is  also  triangular  for  (E,  ft)  e  A  AEH,  with 
vertices  identical  to  those  in  Fig.  11  with  0  — ►  —0. 

We  have  now  described  all  of  Fig.  9  below  the  line  E  =  1/2  (that  is,  below  line  EB).  We 
next  consider  the  situation  for  E  >  1/2.  This  means  that  we  have  to  start  taking  into  account  the 
inequality 

a  <  ^  (129) 

of  Eq.  (126).  This  is  a  constant  upper  bound  on  a  which  will  become  less  than  the  a  coordinate 
of  vertex  (c)  in  Figs.  10  and  11  when  E  >  1/2.  This  results  in  a  horizontal  truncation  of  the 
quadrilateral  of  Fig.  10  so  that  it  becomes  a  pentagon,  and  of  the  triangle  of  Fig.  11  so  that  it 
becomes  a  quadrilateral.  These  are  shown  in  Figs.  12  and  13,  respectively.  The  top  vertices  (/)  and 
( g )  in  both  of  these  figures  are  at  the  upper  bound  a  =  (2  —  E)/ 3,  and  (2  —  e)/3  —  E  =p  (20  —  v)  or 
0  —  (3v±4e=F2)/6,  respectively;  these  vertices  replace  vertex  (c)  in  Figs.  10  and  11.  The  expressions 
for  the  coordinates  of  vertices  (a),  (6),  ( d )  and  (e)  are  identical  to  those  derived  previously.  The 
vertex  coordinates  are  given  in  the  tables  in  their  corresponding  figures. 

Yet  again,  we  ask  for  what  range  of  E  and  n  the  boundaries  pictured  in  Figs.  12  and  13  are 
valid.  As  £  increases,  the  upper  bound  on  a  decreases  until  vertices  (6)  and  (/)  coincide.  This 
happens  when  ( E +fr)/3  =  (2  —  E)/3,  or  £  =  1  —  f/2.  This  is  line  FB  in  Fig.  9.  Above  this  line,  the 
pentagonal  region  of  Fig.  12  degenerates  to  a  quadrilateral,  and  the  quadrilateral  region  of  Fig.  13 
degenerates  to  a  triangle.  In  both  cases,  vertices  ( b )  and  (/)  are  replaced  by  a  new  vertex  (h) 
whose  coordinates  are  a  —  0  —  (2  —  e)/3.  These  situations  are  shown  in  Figs.  14  and  15,  along 
with  corresponding  tables  of  vertex  coordinates. 

The  quadrilateral  of  Fig.  14  will  degenerate  into  the  triangle  of  Fig.  15  when  vertices  (a)  and 
(d)  merge  and  are  replaced  by  vertex  (e).  As  before,  this  degeneracy  happens  when  E  =  n,  and 
this  is  the  boundary  line  CJ  in  the  illustration  of  the  (e,  7f)  plane  shown  in  Fig.  9.  The  triangle  of 
Fig.  15,  in  turn,  degenerates  into  the  empty  set  when  (2  —  E)/3  =  7r  -  e,  or  E  =  3fr/2  -  1.  Referring 
to  Fig.  9,  this  is  the  boundary  line  BC\  thus,  we  see  that  the  quadrilaterals  of  Fig.  14  are  obtained 
when  (e,n)  €  ACFJ,  and  the  triangles  of  Fig.  15  are  obtained  when  (e,  7r)  e  A BCJ. 

Finally,  there  is  one  other  way  that  the  quadrilaterals  of  Fig.  14  can  degenerate.  Vertices  ( d ) 
and  ( g )  will  coincide  if  (2  —  E)/ 3  —  (E  —  fr)/3,  or  E  —  1  +  n/2.  For  values  of  E  greater  than  this, 
vertices  (d)  and  ( g )  are  replaced  by  vertex  (i)  with  coordinates  given  by  a  =  —0  —  (2  —  E)/ 3.  The 
resulting  isoceles  triangular  region  is  shown  in  Fig.  16,  along  with  a  corresponding  table  of  vertex 
coordinates.  This  triangluar  region  degenerates  to  the  empty  set  only  when  E  —  2.  Referring  to 
Fig.  9,  we  see  that  the  triangles  of  Fig.  16  are  obtained  when  (E,  tt)  6  ACFL;  in  fact,  since  these 
regions  are  symmetric  in  0,  they  are  also  symmetric  in  n,  so  their  description  is  the  same  for  all 
(e,  7f)  e  ACDF. 
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Figure  10:  Bounds  on  a  and  /3  for  7r  =  1/5  and  e  —  1/3,  and  the  coordinates  of  the  vertices  given 
as  general  functions  of  e  and  n. 
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Figure  11:  Bounds  on  a  and  /3  for  7r  =  1/3  and  e  =  1/4,  and  the  coordinates  of  the  vertices  given 
as  general  functions  of  e  and  7f. 
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Figure  12:  Bounds  on  a  and  /?  for  n  =  1/5  and 
as  general  functions  of  e  and  #. 


=  2/3,  and  the  coordinates  of  the  vertices  given 
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Figure  13:  Bounds  on  a  and  (3  for  n  =  2/3  and  e  =  11/20,  and  the  coordinates  of  the  vertices 
given  as  general  functions  of  e  and  n. 


51 


BetaBar 


Vertex 

a 

0 

(a) 

0 

0 

(h) 

2  —  £ 

3 

2—e 

3 

(9) 

2— £ 

3 

3vr—4e+2 

6 

(d) 

£  —  7V 

3 

—  £+7f 

3 

Figure  14:  Bounds  on  a  and  /3  for  7r  =  2/3  and  e  =  1,  and  the  coordinates  of  the  vertices  given 
as  general  functions  of  e  and  it. 
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Figure  15:  Bounds  on  a  and  /?  for  n  =  1  and  e  —  3/4,  and  the  coordinates  of  the  vertices  given 
as  general  functions  of  e  and  it. 
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Figure  16:  Bounds  on  a  and  /3  for  n  =  1/2  and  e  =  3/2,  and  the  coordinates  of  the  vertices  given 
as  general  functions  of  e  and  7f. 

B  Basis  Bras  and  Kets  for  FCHC  Model 

In  this  appendix,  we  present  one  possible  choice  for  the  bras  and  kets  of  the  FCHC  model  for 
three-dimensional  fluid  dynamics.  The  basis  bras  are  shown  in  Table  4,  and  the  basis  kets  are 
shown  in  Table  5. 
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Table  4:  Bras  for  the  three-dimensional  projection  of  the  FCHC  model.  The  first  four  rows  are  the 
hydrodynamic  bras,  and  the  last  14  present  one  choice  for  the  kinetic  bras. 


/  -16 

8 

8 

-16 

8 

8 

-16 

8 

8 

-16 

8 

8 

-16 

8 

8 

-16 

8 

8  \ 

12 

0 

-24 

12 

0 

0 

12 

0 

O 

12 

0 

-12 

12 

0 

0 

12 

0 

-12 

-24 

12 

36 

-24 

12 

-12 

-24 

12 

-12 

-24 

12 

12 

-24 

12 

-12 

-24 

12 

12 

-6 

6 

6 

-6 

18 

-18 

-6 

6 

-6 

-6 

0 

6 

-6 

6 

-6 

-6 

0 

6 

-36 

12 

12 

12 

-12 

12 

-36 

12 

12 

-60 

24 

12 

-36 

12 

12 

-60 

24 

12 

0 

0 

0 

-48 

36 

12 

0 

0 

0 

24 

-18 

-6 

0 

0 

0 

24 

-18 

-6 

29 

-7 

-19 

5 

11 

5 

5 

-7 

-1 

29 

-16 

-13 

29 

-7 

-1 

29 

-16 

-13 

-14 

10 

10 

10 

-14 

-2 

10 

-14 

-2 

-14 

10 

4 

-14 

10 

-2 

-14 

10 

4 

-6 

-6 

18 

-6 

6 

-6 

-6 

18 

-18 

-6 

0 

6 

-6 

-6 

6 

-6 

0 

6 

18 

-6 

-6 

-6 

6 

-6 

18 

-6 

6 

18 

-12 

-6 

18 

-6 

-18 

42 

-12 

-6 

-9 

3 

15 

15 

-15 

-9 

-9 

3 

-3 

-9 

0 

9 

-9 

3 

-3 

-33 

24 

9 

36 

-12 

-60 

36 

-12 

12 

36 

-12 

12 

36 

0 

-36 

36 

-12 

12 

36 

-24 

-12 

22 

-2 

22 

-26 

22 

-14 

-2 

-2 

-14 

22 

-14 

16 

-2 

-2 

-14 

22 

-14 

-8 

-28 

20 

20 

20 

-28 

-4 

-4 

-4 

-4 

-28 

20 

8 

-4 

-4 

-4 

-28 

20 

8 

24 

-24 

-24 

24 

0 

0 

24 

0 

0 

24 

-12 

-12 

24 

0 

0 

24 

-12 

-12 

12 

0 

0 

-12 

0 

0 

12 

0 

0 

12 

-12 

0 

12 

0 

0 

12 

-12 

0 

-8 

4 

4 

16 

-8 

-8 

-8 

4 

4 

-8 

10 

-2 

-8 

4 

4 

-8 

10 

-2 

V  12 

-12 

-12 

12 

-12 

12 

12 

-12 

12 

12 

0 

-12 

12 

-12 

12 

12 

0 

-12  / 

Table  5:  Kets  for  the  three-dimensional  projection  of  the  FCHC  model,  corresponding  to  the  choice 
of  bras  in  Table  4.  The  first  four  columns  are  the  hydrodynamic  kets,  and  the  last  14  are  the  kinetic 
kets. 
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14.  ABSTRACT 


We  present  a  general  methodology  for  constructing  lattice  Boltzmann  models  of  hydrody¬ 
namics  with  certain  desired  features  of  statistical  physics  and  kinetic  theory.  We  show  how  a 
methodology  of  linear  programming  theory,  known  as  Fourier- Motzkin  elimination,  provides  an 
important  tool  for  visualizing  the  state  space  of  lattice  Boltzmann  algorithms  that  conserve  a 
given  set  of  moments  of  the  distribution  function.  We  show  how  such  models  can  be  endowed 
with  a  Lyapunov  functional,  analogous  to  Boltzmann’s  H,  resulting  in  unconditional  numeri¬ 
cal  stability.  Using  the  Chapman-Enskog  analysis  and  numerical  simulation,  we  demonstrate 
that  such  entropically  stabilized  lattice  Boltzmann  algorithms,  while  fully  explicit  and  perfectly 
conservative,  may  achieve  remarkably  low  values  for  transport  coefficients,  such  as  viscosity. 
Indeed,  the  lowest  such  attainable  values  are  limited  only  by  considerations  of  accuracy,  rather 
than  stability.  The  method  thus  holds  promise  for  high-Reynolds  number  simulations  of  the 
Navier-Stokes  equations. 
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